%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