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

BHT_operators(r,z,K1,K2,P1,P2,dt,zz)
%% Authored by Joshua Soneson 2007
function[A,B] = BHT_operators(r,z,K1,K2,P1,P2,dt,zz)
% builds operators for IRK2 integration

J = length(r);
M = length(z);

dr = r(2);
dz = z(2);

JM = J*M;
m_ = round(zz*M);

% build matrix operator's vector "bands" 
alpha_plus = ones(JM,1)/dz/dz;
alpha_minus = ones(JM,1)/dz/dz;
bp = zeros(J,1);
bm = zeros(J,1);
bp(1) = 2/dr/dr;
bp(2:J) = (1/dr + 0.5./r(2:J))/dr;
bm(2:J) = (1/dr - 0.5./r(2:J))/dr;
beta_plus = zeros(JM,1);
beta_minus = zeros(JM,1);
for jm=1:JM
  if(mod(jm,M)==0) alpha_plus(jm) = 0; end
  if(mod(jm,M)==1) alpha_minus(jm) = 0; end
  n = ceil(jm/M);
  beta_plus(jm) = bp(n);
  beta_minus(jm) = bm(n);
end
gamma = -2*(1/dr/dr + 1/dz/dz)*ones(JM,1);

% build diffusivity coefficient matrix:
k1 = K1*ones(m_,1);
k2 = K2*ones(M-m_,1);
k = [k1;k2];
bigk = k;
for j=1:J-1
  bigk = [bigk;k];
end
K = spdiags(bigk,0,JM,JM);

% build perfusion coefficeint matrix:
p1 = P1*ones(m_,1);
p2 = P2*ones(M-m_,1);
p = [p1;p2];
bigp = p;
for j=1:J-1
  bigp = [bigp;p];
end
P = spdiags(bigp,0,JM,JM);

% put it all together
A = spdiags([beta_plus alpha_plus gamma alpha_minus beta_minus],[-M -1 0 1 M],JM,JM);
A = A';
A = A*K - P;
B = speye(JM) - 0.5*dt*A;
  

Contact us