Code covered by the BSD License  

Highlights from
High intensity focused ultrasound simulator

High intensity focused ultrasound simulator

by

 

Simulates high intensity focused ultrasound beams and heating effects in layered media

...
%% Authored by Joshua Soneson 2007, updated 2010 %%%%%%%%%%%%%%%%%%%
function[p0,c1,c2,rho1,rho2,N1,N2,G1,G2,gamma1,gamma2,a,b,d,f,R,Z,z_,K] = ...
KZK_parameters(K)
%% Input file for both axisymmetricKZK.m and axisymmetricBHT.m
%% Consists of user-defined input parameter definitions, including 
%% material and transducer parameters, as well as the size of the 
%% computational domain; returns required computed parameters for 
%% integration of the model equations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% User-defined input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% material 1 parameters:
c1 = 1482;		% small-signal sound speed (m/s)
rho1 = 1000;		% mass density (kg/m^3)
alpha1 = .217;		% attenuation at 1MHz (dB/m)
eta1 = 2;		% power of attenuation vs frequency curve
beta1 = 3.5;		% nonlinear parameter
z_ = 5;		% material transition distance (cm)

% material 2 parameters:
c2 = 1629;
rho2 = 1000;
alpha2 = 58;
eta2 = 1;
beta2 = 4.5;

% transducer parameters:
a = 2.5;		% outer radius (cm)
b = 1;			% inner radius (cm)
d = 8;		% focusing depth (cm)
f = 1.5e6;	% frequency (Hz)
P = 100;		% power (W)

% computational domain size:
R = a;			% max radius (cm)
Z = 1.5*d;		% max axial distance (cm)
K = 128;		% number of harmonics included in simulation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computed equation coefficients: %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p0 = sqrt(2*rho1*c1*P/pi/((a/100)^2-(b/100)^2));	% peak pressure at transducer face	
N1 = 2*pi*p0*beta1*(d/100)*f/rho1/c1^3;	% nonlinear coefficient
N2 = 2*pi*p0*beta2*(d/100)*f/rho2/c2^3;
G1 = pi*(a/100)^2*f/c1/(d/100);		% linear pressure gain
G2 = pi*(a/100)^2*f/c2/(d/100);
gamma1 = zeros(K,1);			% attenuation coefficients
gamma2 = zeros(K,1);
h = f*[1:K]'/1e6;
alpha1 = (d/100)*alpha1/8.686;    	  % convert alpha from dB/m to Np
alpha2 = (d/100)*alpha2/8.686;
if(eta1==1)
  gamma1 = alpha1*h.*(1-2*i*log(h)/pi);
elseif(eta1==2)
  gamma1 = alpha1*h.^2;
else
  gamma1 = alpha1*h.^eta1 - 2*i*(alpha1*h.^eta1-alpha1*h)/(eta1-1)/pi;
end
if(eta2==1)
  gamma2 = alpha2*h.*(1-2*i*log(h)/pi);
elseif(eta2==2)
  gamma2 = alpha2*h.^2;
else
  gamma2 = alpha2*h.^eta2 - 2*i*(alpha2*h.^eta2-alpha2*h)/(eta2-1)/pi;
end

% nondimensionalize grid dimensions:
R = R/a;
Z = Z/d;

Contact us