%% 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)