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

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_radial_plots(r,Ir,H(:,round((M+1)/Z)),p5r,p0,rho2,c2,R,a); 
    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_radial_plots(r,Ir,H(:,round((M+1)/Z)),p5r,p0,rho2,c2,R,a); 
    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
KZK_radial_plots(r,Ir,H(:,round((M+1)/Z)),p5r,p0,rho2,c2,R,a); 
KZK_axial_plots(z,Ix,p5x,H(1,:),peak,trough,p0,rho1,rho2,c1,c2,d,Z,a,m_t)

Contact us