%  Calculate the function determining velocity profiles
%  for the common femoral arterie
%
%  Version 1.1, August 25, 1997, JAJ
%    bessel replaced with besselj
%  Version 3.0, 18/9-2017, JAJ
%     Improved version with better display and distances on

%  Physical parameters

rho=1.06e3;                 %  Density of blood   [kg/m^3]
mu=0.004;                   %  Viscosity of blood [kg/m s]
R=4.2/1000;                 %  Radius of vessel   [m]
omega_0=2*pi*62/60;         %  Angular frequency, heart rate = 62 beats/min
V0=0.15;                    %  Mean velocity      [m/s]
deltaT=2*pi/(omega_0*20);   %  Sampling interval of waveform
deltaR=0.05;                %  Radial sampling interval
artery_name='Profiles for Femoral Artery';  %  Name for plots

%   Calculate psi

axes(hp1)
cla
axis([0 1 0 1])
axis off
text (0.3,0.4,['Calculating harmonic components'])
drawnow

psi=0;
for p=1:8
  disp(['Harmonic number ',num2str(p)])
  omega=p*omega_0;
  tau_alpha=sqrt(-1)^(3/2)*R*sqrt(omega*rho/mu);
  index=1;
  for r_rel=0:deltaR:0.99
    Be=tau_alpha*besselj(0,tau_alpha);
    psi(index,p)=(Be-tau_alpha*besselj(0,r_rel*tau_alpha))/(Be-2*besselj(1,tau_alpha));
    index=index+1;
    end
  end


%   Data for phase and amplitude of waveform

Vp=[1.00    0.00
    1.89   32 
    2.49   85 
    1.28  156
    0.32  193 
    0.27  133 
    0.32  155
    0.28  195 
    0.01  310];

%  Convert the data to the correct amplitude

Vp(:,1)=Vp(:,1)*V0;
Vp(:,2)=Vp(:,2)*pi/180;


%  Make the summation for all radial positions and time

profiles=0;
time_index=1;
for time=0:deltaT:2*pi/omega_0
  index=1;
  cla
  axis off
  text (0.3,0.4,['Making flow data. Processed ',num2str(floor(time*62/60*100)),' %'])
  drawnow
  for r_rel=0:deltaR:0.99
    prof=2*V0*(1-r_rel^2);
    for p=1:8
      prof=prof + Vp(p+1,1)*abs(psi(index,p))*cos(p*omega_0*time - Vp(p+1,2) + angle(psi(index,p)));
      end;
    profiles(index,time_index)=prof;
    index=index+1;
    end;
  time_index=time_index+1;
  end;

%  Pad with zeros for the largest radius

profiles(floor(sqrt(2)/deltaR+5.5),time_index-1)=0;

