%% Authored by Joshua Soneson 2007, updated 2010 %%%%%%%%%%%%%%%%%%%
function[p0,c1,c2,rho1,rho2,N1,N2,G1,G2,gamma1,gamma2,a,b,d,f,R,Z,z_,K] = ...
KZK_parameters(K)
%% Input file for both axisymmetricKZK.m and axisymmetricBHT.m
%% Consists of user-defined input parameter definitions, including
%% material and transducer parameters, as well as the size of the
%% computational domain; returns required computed parameters for
%% integration of the model equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% User-defined input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material 1 parameters:
c1 = 1482; % small-signal sound speed (m/s)
rho1 = 1000; % mass density (kg/m^3)
alpha1 = .217; % attenuation at 1MHz (dB/m)
eta1 = 2; % power of attenuation vs frequency curve
beta1 = 3.5; % nonlinear parameter
z_ = 5; % material transition distance (cm)
% material 2 parameters:
c2 = 1629;
rho2 = 1000;
alpha2 = 58;
eta2 = 1;
beta2 = 4.5;
% transducer parameters:
a = 2.5; % outer radius (cm)
b = 1; % inner radius (cm)
d = 8; % focusing depth (cm)
f = 1.5e6; % frequency (Hz)
P = 100; % power (W)
% computational domain size:
R = a; % max radius (cm)
Z = 1.5*d; % max axial distance (cm)
K = 128; % number of harmonics included in simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computed equation coefficients: %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p0 = sqrt(2*rho1*c1*P/pi/((a/100)^2-(b/100)^2)); % peak pressure at transducer face
N1 = 2*pi*p0*beta1*(d/100)*f/rho1/c1^3; % nonlinear coefficient
N2 = 2*pi*p0*beta2*(d/100)*f/rho2/c2^3;
G1 = pi*(a/100)^2*f/c1/(d/100); % linear pressure gain
G2 = pi*(a/100)^2*f/c2/(d/100);
gamma1 = zeros(K,1); % attenuation coefficients
gamma2 = zeros(K,1);
h = f*[1:K]'/1e6;
alpha1 = (d/100)*alpha1/8.686; % convert alpha from dB/m to Np
alpha2 = (d/100)*alpha2/8.686;
if(eta1==1)
gamma1 = alpha1*h.*(1-2*i*log(h)/pi);
elseif(eta1==2)
gamma1 = alpha1*h.^2;
else
gamma1 = alpha1*h.^eta1 - 2*i*(alpha1*h.^eta1-alpha1*h)/(eta1-1)/pi;
end
if(eta2==1)
gamma2 = alpha2*h.*(1-2*i*log(h)/pi);
elseif(eta2==2)
gamma2 = alpha2*h.^2;
else
gamma2 = alpha2*h.^eta2 - 2*i*(alpha2*h.^eta2-alpha2*h)/(eta2-1)/pi;
end
% nondimensionalize grid dimensions:
R = R/a;
Z = Z/d;