k-Wave User Forum » Topic: Differences in focus pressure between k-Wave and analytical transducer model
http://www.k-wave.org/forum/topic/differences-in-focus-pressure-between-k-wave-and-analytical-transducer-model
Support for the k-Wave MATLAB toolboxen-USTue, 18 Jan 2022 01:35:48 +0000http://bbpress.org/?v=1.0.2<![CDATA[Search]]>q
http://www.k-wave.org/forum/search.php
Deepak Sonker on "Differences in focus pressure between k-Wave and analytical transducer model"
http://www.k-wave.org/forum/topic/differences-in-focus-pressure-between-k-wave-and-analytical-transducer-model#post-8283
Fri, 06 Aug 2021 20:03:18 +0000Deepak Sonker8283@http://www.k-wave.org/forum/<p>Thank you I have got the correct pressure with kWaveArray class.
</p>Bradley Treeby on "Differences in focus pressure between k-Wave and analytical transducer model"
http://www.k-wave.org/forum/topic/differences-in-focus-pressure-between-k-wave-and-analytical-transducer-model#post-8254
Wed, 21 Jul 2021 15:13:41 +0000Bradley Treeby8254@http://www.k-wave.org/forum/<p>Hi deepak11,</p>
<p>There are some examples of doing this comparison in the new <a href="http://www.k-wave.org/forum/topic/alpha-version-of-kwavearray-off-grid-sources">kWaveArray</a> class that might show you where you're going wrong. Note, if you're using a grid-based source, you'll have to include an extra scaling factor (or weight the source amplitudes) as discuss in <a href="http://bug.medphys.ucl.ac.uk/papers/2016-Martin-IEEETUFFC.pdf">this paper</a>.</p>
<p>Brad
</p>deepak11 on "Differences in focus pressure between k-Wave and analytical transducer model"
http://www.k-wave.org/forum/topic/differences-in-focus-pressure-between-k-wave-and-analytical-transducer-model#post-8244
Tue, 20 Jul 2021 22:51:29 +0000deepak118244@http://www.k-wave.org/forum/<p>Hi,<br />
First of all thank you for building this amazing toolbox.</p>
<p>I have simulated the bowel transducer in 3D and then compared the focus pressure amplitude with the analytical model of a concave spherical transducer (from Theory of Focusing Radiators by H.T. O'Neil). I also checked this analytical model results from the Sonic Concepts transducer (H-276) and it was correct. I couldn't find the mistake in my code I think I have done everything correctly but the results shouldn't deviate that much (like the peak pressure at the focus must be 58 MHz but the simulated pressure is 46.8 MHz ). I have posted the code below:</p>
<p>(`% define the grid parameters 3D<br />
Nx = 150; % [grid points]<br />
Ny = 150; % [grid points]<br />
Nz = 92;<br />
dx = 0.1e-3; % [m]<br />
dy = 0.1e-3; % [m]<br />
dz = 0.1e-3; % [m]</p>
<p>% create the computational grid<br />
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);</p>
<p>% define the properties of the propagation medium water</p>
<p>medium.sound_speed = 1500; % [m/s]<br />
medium.density = 1000; % [kg/m^3]<br />
medium.alpha_coeff = 0.025; % [dB/(MHz^y cm)]<br />
medium.alpha_power = 1.5;</p>
<p>% define parameters<br />
grid_size = [150,150,92];<br />
bowl_pos = [Nx/2, Ny/2, 2];<br />
diameter = 15e-3; % [m]<br />
radius = 8.5e-3; % [m]<br />
focus_pos = [Nx/2, Ny/2, 87];<br />
freq = 1.5e6; % [Hz]<br />
amp = 2.06e6; % [Pa] </p>
<p>% create bowl<br />
source.p_mask = makeBowl(grid_size, bowl_pos, round(radius / dx), round(diameter / dx)+1, focus_pos, 'Plot', true);</p>
<p>% calculate the time step using an integer number of points per period<br />
ppw = max(medium.sound_speed(:)) / (freq * dx); % points per wavelength<br />
cfl = 0.3; % cfl number<br />
ppp = ceil(ppw / cfl); % points per period<br />
T = 1 / freq; % period [s]<br />
dt = T / ppp; % time step [s]</p>
<p>% calculate the number of time steps to reach steady state<br />
t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2 ) / max(medium.sound_speed(:));<br />
Nt = round(t_end / dt);</p>
<p>% create the time array<br />
kgrid.setTime(Nt, dt);</p>
<p>% define the input signal<br />
source.p = createCWSignals(kgrid.t_array, freq, amp, 0);</p>
<p>% set the sensor mask to cover the entire grid<br />
sensor.mask = ones(Nx, Ny, Nz);<br />
sensor.record = {'p','p_rms','u','u_rms', 'I','I_avg'};</p>
<p>% record the last 3 cycles in steady state<br />
num_periods = 3;<br />
T_points = round(num_periods * T / kgrid.dt);<br />
sensor.record_start_index = Nt - T_points + 1;</p>
<p>% set the input arguements<br />
input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ...<br />
'off', 'PlotScale', [-1, 1] * amp};</p>
<p>% run the acoustic simulation<br />
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});</p>
<p>% convert the absorption coefficient to nepers/m<br />
alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ...<br />
(2 * pi * freq).^medium.alpha_power;</p>
<p>% extract the pressure amplitude at each position<br />
p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq);</p>
<p>% reshape the data, and calculate the volume rate of heat deposition<br />
p = reshape(p, Nx, Ny, Nz);
</p>