Code covered by the BSD License

# High intensity focused ultrasound simulator

### Josh Soneson (view profile)

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

axisymmetricKZK.m
```%% Authored by Joshua Soneson 2008, updated 2010
function[z,r,H,I,Ppos,Pneg]=axisymmetricKZK()
% Driver for axisymmetric KZK integrator.

% get system parameters:
[p0,c1,c2,rho1,rho2,N1,N2,G1,G2,gamma1,gamma2,a,b,d,f,R,Z,z_,K] = KZK_parameters();
K2 = 2*K;

% print coefficients to screen:
fprintf('\n\tp0 = %2.2f MPa\n',1e-6*p0)
fprintf('\tN1 = %1.2f\tN2 = %1.2f\n',N1,N2)
fprintf('\tG1 = %3.2f\tG2 = %3.2f\n\n',G1,G2)

F=0.5*d/a;
if(F<1.37)
fprintf('\tWarning--f/%1.2f exceeds the conditions\n',F)
fprintf('\tunder which KZK is derived (> f/1.37).\n\n')
end

% grid set-up:
[M,J,J_,dz,dr,z,r]=computational_grid(Z,R,max(G1,G2),a,d,gamma2(1),N2);

% dependent variables:
u = zeros(2*J,K);
w = zeros(2*J,K);
limit = 1/sqrt(1-(a/d)^2);
v = initial_condition(J,K,G1,r,b*limit/a,limit);
v(1:J,1) = v(1:J,1).*sqrt(1-(r/d).^2);
v(J+1:2*J,1) = v(J+1:2*J,1).*sqrt(1-(r/d).^2);

% set up discretization operators:
for k=1:K
[A1(k).IRK1,A1(k).IRK2,A1(k).CN1,A1(k).CN2] ...
= KZK_operators(r,R,G1,dz,dr,J,k,gamma1(k));
[A2(k).IRK1,A2(k).IRK2,A2(k).CN1,A2(k).CN2] ...
= KZK_operators(r,R,G2,dz,dr,J,k,gamma2(k));
end
k1 = zeros(2*J,1);	% IRK slope vectors
k2 = zeros(2*J,1);
b1 = 1/sqrt(2);		% IRK coefficients
b2 = 1 - b1;

% parameters for nonlinear integration:
mu1 = N1*K*dz/pi;		% nonlinear term integration parameters
mu2 = N2*K*dz/pi;
cutoff1 = gamma1(1)/10/N1;	% cutoffs for performing nonlinear integration
cutoff2 = gamma2(1)/10/N2;
X = zeros(1,K2);		% data vectors
Y = zeros(1,K2);
Xpeak = zeros(1,K2);

% for plotting routines:
H = zeros(J_,M);		% Heating rate matrix
H2 = zeros(J_,M);
H(:,1) = real(gamma1(1))*(v(1:J_,1).^2 + v(J+1:J+J_,1).^2);
I = zeros(J_,M);
I(:,1) = v(1:J_,1).^2 + v(J+1:J+J_,1).^2;
I_td = zeros(J,1);
dt = 1/f/(K2-1);
Ix = zeros(M,1);		% axial intensity vector
Ix(1) = v(1,1)^2 + v(J+1,1)^2;
Ir = zeros(J_,1);		% radial-focal intensity vector
if(K<5) kk = K; else kk = 5; end
p5x = zeros(kk,M); 		% first (up to) 5 harmonic pressure amplitudes
p5x(1,1) = sqrt(v(1,1)*v(1,1) + v(J+1,1)*v(J+1,1));
p5r = zeros(kk,J_);
peak = zeros(M,1);		% waveform peak
trough = zeros(M,1);		% waveform trough
peak(1) = 1;
trough(1) = -1;

% for monitoring simulation:
amp_K = 0;			% amplitude of Kth harmonic

% for determining peak pressure waveform and location of its occurence:
m_t = ceil(z_/dz/d);		% index of last meshpoint in material 1
if(m_t>M) m_t = M; end
m_f = round(M/Z);		% index of meshpoint nearest focus
p_peak = 0;			% peak axial pressure
z_peak = 0;			% distance at which peak axial pressure occurs
Ppos = zeros(J,M);
Pneg = zeros(J,M);

% integrate the equations through material 1:
fprintf('\tIntegrating KZK equation...\n')
fprintf('\t\tz (cm)\ttime (hr:min:sec)\n')
t_start = clock;
p1 = 0;
for m=1:m_t-1
w = v;
[v,X,Ppos(:,m),Pneg(:,m),I_td] = ...
TDNL(v,X,Y,K,J,mu1,cutoff1,Ppos(:,m),Pneg(:,m),I_td);
I_td = f*dt*I_td/dz;
peak(m+1) = max(X);
trough(m+1) = min(X);
if(m==m_f)
if(K>1)
Ir = sum((v(1:J_,:)').^2);
Ir = Ir + sum((v(J+1:J+J_,:)').^2);
else
Ir = v(1:J_).^2 + v(J+1:J+J_).^2;
end
p5r = sqrt(v(1:J_,1:kk).^2 + v(J+1:J+J_,1:kk).^2);
end
if(peak(m+1)>p_peak)
p_peak = peak(m+1);
z_peak = z(m+1);
Xpeak = X;
end
if(z(m)<0.3)
for k=1:K
k1 = A1(k).IRK1 \ (A1(k).IRK2*v(:,k));
k2 = A1(k).IRK1 \ (A1(k).IRK2*(v(:,k) + dz*b1*k1));
u(:,k) = v(:,k) + dz*(b1*k1 + b2*k2);
end
else
for k=1:K
u(:,k) = A1(k).CN1 \ (A1(k).CN2*v(:,k));
end
end
v = u;
pl = sqrt(v(1,K).^2 + v(J+1,K).^2);
if(pl > amp_K) amp_K = pl; end
for j=1:J_
H(j,m+1) = sum(real(gamma1(:)').*(u(j,:).^2 + u(J+j,:).^2));
H2(j,m+1) = I_td(j);
I(j,m+1) = sum(u(j,:).^2 + u(J+j,:).^2);
end
Ix(m+1) = sum(v(1,:).^2 + v(J+1,:).^2);
if(Ix(m+1)>2*G1*G1)
fprintf('\tStopped - computation became unstable at z = %2.1f cm.\n',d*z(m))
r = a*r(1:J_);
z = d*z;
KZK_axial_plots(z,Ix,p5x,H(1,:),peak,trough,p0,rho1,rho2,c1,c2,d,Z,a,m_t)
return
end
p5x(:,m+1) = sqrt(v(1,1:kk).*v(1,1:kk)+v(J+1,1:kk).*v(J+1,1:kk));
p2 = floor(10*(m+1)/M);
p1 = timing(p1,p2,t_start,z(m+1),d);
end

% material 2:
v = 2*rho2*c2*v/(rho1*c1+rho2*c2);
for m=m_t:M-1
w = v;
[v,X,Ppos(:,m),Pneg(:,m),I_td] = ...
TDNL(v,X,Y,K,J,mu2,cutoff2,Ppos(:,m),Pneg(:,m),I_td);
I_td = f*dt*I_td/dz;
peak(m+1) = max(X);
trough(m+1) = min(X);
if(m==m_f)
if(K>1)
Ir = sum((v(1:J_,:)').^2);
Ir = Ir + sum((v(J+1:J+J_,:)').^2);
else
Ir = v(1:J_).^2 + v(J+1:J+J_).^2;
end
p5r = sqrt(v(1:J_,1:kk).^2 + v(J+1:J+J_,1:kk).^2);
end
if(peak(m+1)>p_peak)
p_peak = peak(m+1);
z_peak = z(m+1);
Xpeak = X;
end
if(z(m)<0.3)
if(m==m_t)
for k=1:K
k1 = A1(k).IRK1 \ (A1(k).IRK2*v(:,k));
k2 = A2(k).IRK1 \ (A2(k).IRK2*(v(:,k) + dz*b1*k1));
u(:,k) = v(:,k) + dz*(b1*k1 + b2*k2);
end
u = (1-(rho1*c1-rho2*c2)/(rho1*c1+rho2*c2))*u;
else
for k=1:K
k1 = A2(k).IRK1 \ (A2(k).IRK2*v(:,k));
k2 = A2(k).IRK1 \ (A2(k).IRK2*(v(:,k) + dz*b1*k1));
u(:,k) = v(:,k) + dz*(b1*k1 + b2*k2);
end
end
else
for k=1:K
u(:,k) = A2(k).CN1 \ (A2(k).CN2*v(:,k));
end
end
v = u;
pl = sqrt(v(1,K).^2 + v(J+1,K).^2);
if(pl > amp_K) amp_K = pl; end
for j=1:J_
H(j,m+1) = sum(real(gamma2(:)').*(u(j,:).^2 + u(J+j,:).^2));
H2(j,m+1) = I_td(j);
I(j,m+1) = sum(u(j,:).^2 + u(J+j,:).^2);
end
Ix(m+1) = sum(v(1,:).^2 + v(J+1,:).^2);
if(Ix(m+1)>2*G2*G2)
fprintf('\tStopped - computation became unstable at z = %2.1f cm.\n',d*z(m))
r = a*r(1:J_);
z = d*z;
KZK_axial_plots(z,Ix,p5x,H(1,:),peak,trough,p0,rho1,rho2,c1,c2,d,Z,a,m_t)
return
end
p5x(:,m+1) = sqrt(v(1,1:kk).*v(1,1:kk)+v(J+1,1:kk).*v(J+1,1:kk));
p2 = floor(10*(m+1)/M);
p1 = timing(p1,p2,t_start,z(m+1),d);
end
H(:,1:m_t) = 1e-4*p0*p0*H(:,1:m_t)/rho1/c1/d;	% dimensionalize H
H(:,m_t+1:M) = 1e-4*p0*p0*H(:,m_t+1:M)/rho2/c2/d;
H2(:,1:m_t) = 1e-4*0.5*p0*p0*H2(:,1:m_t)/rho1/c1/d;	% dimensionalize H
H2(:,m_t+1:M) = 1e-4*0.5*p0*p0*H2(:,m_t+1:M)/rho2/c2/d;
H = real(H + H2);
I(:,1:m_t) = 1e-4*0.5*p0*p0*I(:,1:m_t)/rho1/c1;	% dimensionalize I
I(:,m_t+1:M) = 1e-4*0.5*p0*p0*I(:,m_t+1:M)/rho2/c2;

fprintf('\tmax(|p_K|/p0) = %e\n',amp_K)

r = a*r(1:J_);	% rescale r and chop so that PML region is excluded in plots
z = d*z;	% rescale z
Ppos = Ppos(1:J_,:);
Pneg = Pneg(1:J_,:);

% call plotting routines:
if(K>1)
plot_waveform(p0,f,Xpeak,d*z_peak,K2);
end