clear all
clc

%% *** Read Pressure Gradient and Perform Fourier Analysis ************
%%load experimental data 
load('Pressure_Gradient_CC.mat')

t = data(:,1); %time 
pz = data(:,2);
Qexp = data(:,3);

period=max(t);
np=length(pz);

figure
plot(-pz)
xlabel('frequency')
ylabel('pz')
frequency=1/period;

pfz=fft(pz)./sqrt(np);
pfz_mod=abs(pfz);
freq=(((1:np)-1)/period)';

figure
scatter(freq,abs(pfz),'.')
xlabel('frequency')
ylabel('fft(pz)')
%%Select number of harmonics 

nharm=21;

pfz_selected=pfz(1:nharm);

%% *** Rheological Parameters and Geometry ************** 
rho = 1060;
r0 = 0.0055;
mu = 0.004;


for i=1:nharm
    omega(i,1)=2*pi*frequency*i;   %pulsation
end
alpha=sqrt(r0^2*omega*rho/mu);


%% *** Calculate the velocity profile using Womersley's theory ************** 
% Get the velocity profile for each harmonic vharm

radial_elts=25;


velocity=zeros(np,radial_elts);

r_mat=linspace(0,r0,radial_elts);    % From m to cm
alpha_upset=alpha(1:length(alpha)-1);
omega_upset=omega(1:length(omega)-1);

for ti=1:np
    c_time=t(ti);
    for ind_r=1:radial_elts   % RADIAL ***********************************************
        r=r_mat(ind_r);
        v0=pfz_selected(1)/(4*mu)*(r0^2-r^2);  % Poiseuille component
        vh=0;
        for m=1:nharm-1
            vharm = real( pfz_selected(m+1)/(1j*rho*omega_upset(m))*(1- besselj(0,(1j^(3/2)*alpha_upset(m)*r/r0))/besselj(0,(1j^(3/2)*alpha_upset(m))))* exp(1j*omega_upset(m)*c_time) );
            vh=vh+vharm;
        end             
        w=v0+vh;   
        
        velocity(ti,ind_r)=w/100;           % m/s
        
    end     % RADIAL
end     % TIME

figure
plot(r_mat,velocity(np/4,:))
xlabel('radius')
ylabel('velocity')

%% *** Calculate flow rate ************** 
%integrate numrically 