hi every one.i want to obtain the rms pressure of a HIFU with continuous excitation in the plane perpendicular to HIFU but i didnt get answer.the medium which i simulated consist of fat ,muscle,rib and liver.here is my code.i appreciate any help.i am looking forward getting answer from you as soon as possible.thanks in advance.
clear all;
h=50;
% simulation settings
DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations
MASK_PLANE = 'yz'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 10; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]
% set total number of grid points not including the PML
Nx = 170- 2*PML_X_SIZE; % [grid points]
Ny = 170- 2*PML_Y_SIZE; % [grid points]
Nz = 170- 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 15e-2; % [m]
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = dx; % [m]
dz = dx; % [m]
% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
medium.sound_speed = 1578*ones(Nx,Ny,Nz); % [m/s]
medium.sound_speed (:,:,1:h)=1430;
medium.sound_speed (:,:,h:h+10)=1430;
medium.sound_speed (:,:,h+10:h+30)=1580;
medium.sound_speed (:,:,h+20:h+30)=1580;
medium.sound_speed (:,1:10,h+20:h+30)=3198;
medium.sound_speed (:,20:30,h+20:h+30)=3198;
medium.sound_speed (:,40:50,h+20:h+30)=3198;
medium.sound_speed (:,60:70,h+20:h+30)=3198;
medium.sound_speed (:,80:90,h+20:h+30)=3198;
medium.sound_speed (:,100:110,h+20:h+30)=3198;
medium.sound_speed (:,120:130,h+20:h+30)=3198;
medium.sound_speed (:,140:150,h+20:h+30)=3198;
medium.density = 1050*ones(Nx,Ny,Nz); % [kg/m^3]
medium.density(:,:,1:h)=928;
medium.density(:,:,h:h+10)=928;
medium.density(:,:,h+10:h+30)=1041;
medium.density(:,1:10,h+20:h+30)=1990;
medium.density(:,20:30,h+20:h+30)=1990;
medium.density(:,40:50,h+20:h+30)=1990;
medium.density(:,60:70,h+20:h+30)=1990;
medium.density(:,80:90,h+20:h+30)=1990;
medium.density(:,100:110,h+20:h+30)=1990;
medium.density(:,120:130,h+20:h+30)=1990;
medium.density(:,140:150,h+20:h+30)=1990;
medium.alpha_coeff = 0.45*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)]
medium.alpha_coeff (:,:,1:h)=.6;
medium.alpha_coeff (:,:,h:h+10)=.6;
medium.alpha_coeff (:,:,h+10:h+30)=.6;
medium.alpha_coeff (:,1:10,h+20:h+30)=3.54;
medium.alpha_coeff (:,20:30,h+20:h+30)=3.54;
medium.alpha_coeff (:,40:50,h+20:h+30)=3.54;
medium.alpha_coeff (:,60:70,h+20:h+30)=3.54;
medium.alpha_coeff (:,80:90,h+20:h+30)=3.54;
medium.alpha_coeff (:,100:110,h+20:h+30)=3.54;
medium.alpha_coeff (:,120:130,h+20:h+30)=3.54;
medium.alpha_coeff (:,140:150,h+20:h+30)=3.54;
medium.alpha_power = 1.5;
medium.BonA =0;
dt_stability_limit = checkStability(kgrid, medium);
% create the time array
t_end =15e-6; % [s]
kgrid.t_array = makeTime(kgrid,medium.sound_speed,.1, t_end);
% DEFINE THE ULTRASOUND TRANSDUCER
% define a curved transducer element
b=spheregen(h);
source.p_mask=b;
% =========================================================================
% define a time varying sinusoidal source
source_freq = 1e6; % [Hz]
source_mag = 0.5; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
% filter the source to remove any high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
% =========================================================================
% DEFINE SENSOR MASK
% =========================================================================
% define a sensor mask through the central plane
sensor.mask = zeros(Nx, Ny, Nz);
switch MASK_PLANE
case 'yz'
% define mask
sensor.mask(Nx/2, :,:) = 1;
% store y axis properties
Nj = Ny;
j_vec = kgrid.y_vec;
j_label = 'y';
case 'xz'
% define mask
sensor.mask(:, Ny/2, :) = 1;
% store z axis properties
Nj = Nz;
j_vec = kgrid.z_vec;
j_label = 'z';
end
% set the record mode such that only the rms and peak values are stored
if USE_STATISTICS
sensor.record = {'p_rms'};
end
voxelPlot(sensor.mask);
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
input_args = {'PlotLayout',true,'PMLInside', false, 'PlotPML',true, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotScale', [-source_mag/2,source_mag/2]};
% stream the data to disk in blocks of 100 if storing the complete time
% history
if ~USE_STATISTICS
input_args = [input_args {'StreamToDisk', 100}];
end
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium,source, sensor, input_args{:});
% voxelPlot(double(source.p_mask | cart2grid(kgrid, sensor.mask)));
% =========================================================================
% COMPUTE THE BEAM PATTERN USING SIMULATION STATISTICS
% =========================================================================
if USE_STATISTICS
% reshape the returned rms and max fields to their original position
sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]);
% plot the beam pattern using the pressure rms
figure;
imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, sensor_data.p_rms/1e6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using RMS Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;
% end the example
return
end
%%below function creates HIFU with different radius
function [a] = sphreregen(h)
a=makeSphere(150,150,150,50);
b=zeros(size(a));
for i=1:101
b(:,:,i)=a(:,:,24+i);
end
b(:,:,h:end)=0;
voxelPlot(b);
end