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

KZK_operators.m
```%% Authored by Joshua Soneson 2007
function[IRK1,IRK2,CN1,CN2]=make_matrices(r,R,G,dz,dr,J,k,gamma)
%% Computes second-order diagonally implicit Runge-Kutta and
%% Crank-Nicolson operators for integrating the kth harmonic
%% equation.

% PML vector
v = zeros(J,1);
for j=1:J
if(r(j)>R) v(j) = 16*(r(j)-R).^2; end
end
v = exp(i*pi*v/2);
e = ones(J,1);

% build finite-difference matrix:
r(1) = 1;       % prevent divide by zero error
supdiag = (v/dr + 0.5./r)/dr;
diagonal = -2*v/dr/dr;
subdiag = (v/dr - 0.5./r)/dr;
A = spdiags([supdiag diagonal subdiag],-1:1,J,J);
A = conj(A');
A(1,2) = -A(1,1);
A = 0.25*A/G/k;

% split into real and imaginary parts:
ReA = real(A);
ImA = imag(A);
alpha = real(gamma)*speye(J);
beta = imag(gamma)*speye(J);

% build real-valued block system:
M = [ImA-alpha,-ReA-beta;
ReA+beta,ImA-alpha];
I = speye(2*J);

% DIRK2 matrices:
IRK1 = I - dz*(1-1/sqrt(2))*M;
IRK2 = M;

% CN matrices:
CN1 = I - 0.5*dz*M;
CN2 = I + 0.5*dz*M;
```

Contact us