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_parameters(z,r,H)
%% Authored by Joshua Soneson 2007
function[C1,C2,k1,k2,w1,w2,N,t,input,dt,J,M,z,r,H,T0] = BHT_parameters(z,r,H)
% determines parameters used in BHT equation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% User-defined input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% material 1:
C1 = 4180;	% heat capacity (J/kg/K)
k1 = 0.6;	% thermal conductivity (W/m/K)
w1 = 0;		% perfusion rate (kg/m^3/K)

% material 2:
C2 = 4180;
k2 = 0.6;
w2 = 20;


% ambient temperature:
T0 = 37;	% (degrees C)

% continuous sonication
t_i = 0.3;	% initial sonication duration (s)

% periodic sonications
n_c = 5;		% number of pulse cycles
D = 20;		% duty cycle (%)
t_p = 0.5;	% pulse cycle period (s)

% cooling period
t_c = 5.2;	% cool-off duration (s)

% grid coarsening
r_skip = 2;	% reduction factor in r-direction
z_skip = 4;	% reduction factor in z-direction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Computed parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% determine proper timestep:
if(n_c==0)
  dt = min([0.1,t_i/5]);
elseif(t_i==0)
  dt = min(0.1,0.01*D*t_p/5);
else
  dt = min([0.1,0.01*D*t_p/5,t_i/5]);
end

% determine number of integration steps at each stage:
N_i = round(t_i/dt);
if(n_c==0)
  N_p=0;
else
  N_p = round(t_p/dt);
end
N_c = round(t_c/dt);
N = N_i+n_c*N_p+N_c;		% total number of integration steps
T = dt*(N-1);			% total simulation duration

% build input vector, which catalogs the HIFU beam on/off operation 
if(N_p~=0)
  pulse = zeros(1,N_p);
  for n=1:round(0.01*D*N_p);
    pulse(n) = 1;
  end
  pulses = pulse;
  for m=2:n_c
    pulses = [pulses,pulse];
  end
end
if(N_c~=0)
  cooloff = zeros(1,N_c);
end 
if(N_i~=0)
  initial = ones(1,N_i);
end
for n=1:N_i
  input(n) = initial(n);
end
for n=N_i+1:N_i+n_c*N_p
  input(n) = pulses(n-N_i);
end
for n=N_i+n_c*N_p+1:N
    input(n) = cooloff(n-N_i-n_c*N_p);
end

% build grid vectors:
t = [0:dt:T];	% time nodes
[J,M] = size(H);
r = r(1:r_skip:round(J/2));
z = z(1:z_skip:M);
H = H(1:r_skip:round(J/2),1:z_skip:M);
[J,M] = size(H);
fprintf('\n\tdt = %2.2f sec\tN = %d\n',dt,N)
fprintf('\tdr = %2.2f mm\tJ = %d\n',10*r(2),J)
fprintf('\tdz = %2.2f mm\tM = %d\n\n',10*z(2),M)

Contact us