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

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;

```