%The UIV code gives the correlation averaged displacement field at
%successive phases throughout the pulse cycle. This function then
%calculates the velocity fields from those displacements. This includes a
%speed=distance/time calculation and correction due to the accelerating
%motion of the scatters within the sweeping ultrasound beam.

%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 [VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, int, fluid,period,split)
%for normal PIV delta=lines (the number of lines used in the whole image)
%for interlaced PIV delta=dleta (the number of lines between the images)
%int=0 normal PIV
%int=1 interlaced PIV
%Vx, Vy are DISPLACEMENTS between the image pairs - they are not
%velocities!!
 
deltaT=(delta)*(2*depth/1540+23.488e-6); %in s
if int==1
    timestep=(2*128+(delta-1))*(2*depth/1540+23.488e-6);
    totalTime=timestep*period;
    dt=totalTime./split;
else
    dt=deltaT;
end
xpix=0.001*38/128;
VxM=Vx.*xpix; %convert x displacement from pix to m
VxM=VxM./deltaT; %convert m displacement to velocity
if fluid=='water'
    ypix=3.7425e-5/2;
elseif fluid=='glyce'
    ypix=2.3653e-05
end
VyM=Vy.*ypix; %convert y displacement from pix to m
VyM=VyM./deltaT; %convert m displacement to velocity

VxC=zeros(size(VxM));
VyC=zeros(size(VyM));

%sweep speed correction
Vs=(38*1e-3/128)/(2*depth/1540+23.488*1e-6); %sweep speed in m/s
if int==1
    Vs=Vs/2; %interlaced
end
VS=ones(size(Vx))*Vs;
[a,b,c]=size(VxM);

for i=1:a
    for j=1:b
        vx=VxM(i,j,:);
        vx=reshape(vx,[],1);
        if (i==30 & j==120)
            test=1;
        end
        vy=VyM(i,j,:);
        vy=reshape(vy,[],1);
        
        En=i*xpix./(12*dt*Vs)*(-1 -vx./Vs - i*xpix./(2*dt*Vs));
        Dn=2*i*xpix./(3*dt*Vs)*(1+vx./Vs + i*xpix./(dt*Vs));
        Cn=1+vx./Vs -5.*(i*xpix./(2*dt*Vs)).^2;
        Bn=2*i*xpix./(3*dt*Vs)*(-1-vx./Vs +i*xpix./(dt*Vs));
        An=i*xpix./(12*dt*Vs)*(1 +vx./Vs - i*xpix./(2*dt*Vs));
        
        K=length(vx);
        M=zeros(K);
        M(1,1)=Cn(1);
        M(1,2)=Dn(1);
        M(1,3)=En(1);
        M(1,K)=Bn(1);
        M(1,K-1)=An(1);
        
        M(2,1)=Bn(1);
        M(2,2)=Cn(1);
        M(2,3)=Dn(1);
        M(2,4)=En(1);
        M(2,K)=An(1);
        
        for k=3:K-2
            M(k,k-2)=An(k);
            M(k,k-1)=Bn(k);
            M(k,k)=Cn(k);
            M(k,k+1)=Dn(k);
            M(k,k+2)=En(k);
        end

        M(K-1,K)=Dn(K);
        M(K-1,K-1)=Cn(K);
        M(K-1,K-2)=Bn(K);
        M(K-1,K-3)=An(K);
        M(K-1,1)=En(K);

        M(K,K)=Cn(K);
        M(K,K-1)=Bn(K);
        M(K,K-2)=An(K);
        M(K,1)=Dn(K);
        M(K,2)=En(K);
        
        vxc2=M\vx;
        VxC(i,j,:)=vxc2;
        vyc2=M\vy;
        VyC(i,j,:)=vyc2;
    end
end