%This function calculates a fully developed velocity profile for a straight
%tube according to Womersley's theory. It is based on knowing the mean
%velocity.
%For full details see "Ultrasound Imaging Velocimetry with interleaved
%images for improved pulsatile arterial flow measurements: a new correction
%method, experimental and in vivo validation", J Royal Soc Interface, 2017
%K H Fraser, C Poelma 2016
function vel= womer_Vmean(t, freqs, V_mean, R, nu)
%t = time
%freqs = number of harmonics to use
%V_peak = peak velocity waveform (same number of points as t)
%R = radius of artery
[p1,q1]=size(t);
tot_size=p1*q1;
q=length(t);
p=tot_size/q;
%period time is just the last time in the time vector
T=t(end);
% main execution loop. calls the function d2 which sets up graphs and calls
% function d1 which actually calculates and sums Womersley profile harmonics.
n=1;
for r=1:q
velprof = d2(t(r),freqs, V_mean, R, T, nu);
vel(r,:) = velprof; %whole profile
end
function velprof = d2(t,freqs, V_mean, R, T,nu)
% the fourier coefficients. Dividing by the number of samples is necessary
% here as the fourier coefficients are calculated with a factor of N.
B=fft(V_mean)/length(V_mean);
% sets up a complex or real grid, depending on what you are plotting
%z = cplxgrid(20);
z=0:0.01:1;
velprof = 2*B(1)*(1-abs(z).^2) + d1(z,t,B,R, freqs, T,nu);
% this function actually calculates the Womersley harmonics. They are summed
% in the variable k
function k = d1(z,t,B,R, freqs, T,nu)
% frequency -in rad/s
w=2*pi/T;
% Womersley number for mode 0
A0=R*sqrt(w/nu);
k=0;
% this sums the various components to form the complete Womersley profile
for y =2:freqs
An=A0*sqrt(y-1);
top = besselj(0,An*1i^(3/2))-besselj(0,An*abs(z)*1i^(3/2)); %
bottom = besselj(0,An*1i^(3/2)) - 2*besselj(1,An*1i^(3/2))/(An*i^(3/2));% in these 2 lines Holdsworth paper has i^-(3/2) but now fairly sure it should be +
k = k+2*B(y)*(top/bottom)*exp(w*i*(y-1)*t);
end