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

axisymmetricBHT(z,r,H)
%% Authored by Joshua Soneson 2007, updated 2010
function[Tmax_mat,Dmat] = axisymmetricBHT(z,r,H)
% Bioheat transfer equation integrator
 
% get parameters:
[p0,c1,c2,rho1,rho2,N1,N2,G1,G2,gamma1,gamma2,a,b,d,f,R,Z,z_,K] = ...
KZK_parameters();

% get more parameters:
[C1,C2,k1,k2,w1,w2,N,t,input,dt,J,M,z,r,H,T0] = BHT_parameters(z,r,H);
D1 = 1e4*k1/C1/rho1;			% diffusivity of material 1
D2 = 1e4*k2/C2/rho2;			% diffusivity of material 2
P1 = w1/rho1;
P2 = w2/rho2;

% create Crank-Nicolson operators for BHT:
[A,B] = BHT_operators(r,z,D1,D2,P1,P2,dt,z_/d/Z);
m_ = round(z_*M/d/Z);			% index @ material interface 
if(m_>M) m_ = M; end
H(:,1:m_) = 1e6*H(:,1:m_)/C1/rho1;	% rescale H to degrees/second
H(:,m_+1:M) = 1e6*H(:,m_+1:M)/C2/rho2;

JM = J*M;
Q = zeros(JM,1);			% heating source vector
T = zeros(JM,1);			% temperature vector
D = zeros(JM,1);			% thermal dose vector
Tmat = zeros(J,M);			% temp matrix
Dmat = zeros(J,M);			% dose matrix
Q = vektorize(Q,H,J,M);			% column-stack the heat matrix
s1 = zeros(JM,1);			% slopes for IRK integrator
s2 = zeros(JM,1);
Tpeak = zeros(N,1);			% temp vs time vector
Tmax_vec = zeros(JM,1);			% max temp distribution vector
Tmax_mat = zeros(J,M);			% max temp distribution matrix
tt = 0;

% Integrate BHT:
fprintf('\tIntegrating BHT equation...\n')
fprintf('\t\tt (sec)\ttime (hr:min:sec)\n')
t_start = clock;
p1 = 0;
for n=1:N-1
  s1 = A*T + input(n)*Q;
  s2 = B \ (A*(T+0.5*dt*s1) + input(n+1)*Q);
  T = T + 0.5*dt*(s1+s2);
  D = D + equivalent_time(T,JM,T0);
  Tpeak(n+1) = max(T);
  if(Tpeak(n+1)>tt)
    tt = Tpeak(n+1);		
    Tmax_vec = T;
    t_peak = t(n+1);
  end
  p2 = floor(10*(n+1)/N);
  p1 = timing(p1,p2,t_start,t(n+1),1);
end

Tmax_mat = matrixize(Tmax_vec,Tmax_mat,J,M);
Dmat = dt*matrixize(D,Dmat,J,M)/60;
BHT_plots(z,r,t,t_peak,Tpeak,Tmax_mat,Dmat,T0)




Contact us