%   Propagate a pulse to simulate the non-linear distortion
%
%  Calling: pulse = propagate_pulse (P, z_max, fs)
%
%  Input:  P     - Pulse to propagate
%          z_max - Distance to propagate
%          fs    -  Sampling frequency
%
%   Version 1.1, April 22, 2003, JAJ

function  pulse = propagate_pulse (P, z_max, fs)

%  Set various parameters for propagation

c=1540;              %  Speed of sound     [m/s]
dz=0.5/1000;         %  Step size          [m]
dt=dz/c;             %  Step in seconds    [s]
K=round(z_max/dz);   %  Number of iterations
rho=1000;            %  Density            [kg/m^3]
Z=rho*c;             %  Characteristic acoustic impedance
BdivA=5.2;           %  B/A parameter

%  Make the propagation of the pulse

n=1:length(P);
for k=1:K
    U=P/Z;
    c_diff = c*(1+BdivA/2*U/c).^(2/BdivA+1) - c;
    shift = dt*c_diff/c*fs;
    P=(ceil(shift)-shift).*P(n+floor(shift)) + (shift-floor(shift)).*P(n+ceil(shift));
end;

pulse=P;
