%This script generates the images for figure 4 and performs associated
%analysis.

%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

for delta=[15 31 47 63 79 95]
%% -----set up------
folder=['../data/flow_phantom/a/UIV/2cm/'];
name=['/phantom_a_128lines_' int2str(delta) 'delta_20000depth7000000freq_out']
load([folder int2str(delta) name '.mat']);
depth=0.02;
ypix=2.3653e-05;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'glyce',period, split);
if delta==7
    indices=[1:50];
    centre=[1.9,1.1];
elseif delta==15
    indices=[47:50 1:46];
    centre=[1.9,1.05];
elseif delta==31
    indices=[4:50 1:3];
    centre=[1.9,1.05];
elseif delta==47
     indices=[25:50 1:24];
    centre=[1.9,1.05];
elseif delta==63
    indices=[14:50 1:13];
    centre=[1.9,1.05];
elseif delta==79
    indices=[26:50 1:25];
    centre=[1.9,1.05];
elseif delta==95
    indices=[19:50 1:18];
    centre=[1.9,1.05];
end

%--sort out the time steps--
depth=0.02;
lines=128;
timestep=(lines*2+delta)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% -----velocity profile---------
clear vprofile
yfrom=100; yto=180;
for t=indices
    vprofile(:,t)=mean(-VxC([1:22 36:63],yfrom:yto,t),1);
end
r=([1:length(vprofile(:,t))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;
%% ----flow rate-comparison----
% ----flow rate-transonics----
M=dlmread('../data/flow_phantom/a/flow_rate/phantom_a_flow.csv',',',[4 0 936 14]);
timeT=M(:,1);
timeT=timeT./1000;
flowT=M(:,12);
flowTMax=M(:,14);
flowTMin=M(:,15);

%scale by bucket flow rate to eliminate speed of sound error with
%Transonics
scaleFactor=383.8./mean(flowT);
flowTMax=flowTMax.*scaleFactor;
flowTMin=flowTMin.*scaleFactor;
flowT=flowT.*scaleFactor;

%---flow rate UIV---
[flow2, Flow, upperEdge, lowerEdge]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 390, 440, 670, 730);
flow=flow2(indices);
figure(100); hold on;
if delta==7
    plot(times,flow,'Color',[1 0 1],'LineWidth',0.1);
elseif delta==15
    plot(times,flow,'m','LineWidth',0.1);
elseif delta==31
    plot(times,flow,'r','LineWidth',0.1);
elseif delta == 47
    plot(times,flow,'Color',[1 0.5 0],'LineWidth',0.1);
elseif delta==63
    plot(times,flow,'g','LineWidth',0.1);
elseif delta==79
    plot(times,flow,'c','LineWidth',0.1);
elseif delta==95
    plot(times,flow,'Color',[0.55 0 0.55],'LineWidth',0.1);
end
Tindices=[30:length(flowT) 1:29];
flowT2 = interp1([timeT], [flowT(Tindices)],times(1:end-1));
meanFlowerror=mean(abs(flow(1:end-1)-flowT2'))
stdFlowerror=std(abs(flow(1:end-1)-flowT2'))
percentFlowerror=meanFlowerror./(max(flow)).*100
stdpercentFlowerror=stdFlowerror./(max(flow)).*100
%% ----womersley profile----
figure;
for t=1:50
    R(t)=ypix*mean(lowerEdge(:,t)-upperEdge(:,t))/2;
end
Rmean=mean(R);
flow_f=flowT(Tindices);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
nu=2.7574e-06;

vmean=flowT(Tindices)./(pi*Rmean^2);
vmean2=vmean./(60*1000*1000);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
vel3= womer_Vmean(timeT(pos), 20, vmean2(pos), Rmean,nu);

vprofile_average=(vprofile(round(length(vprofile)/2):end,:)+vprofile(fliplr(1:round(length(vprofile)/2)),:))/2;

figure;  plot([1:50], real(vel3(:,1)),'k');
hold on; plot([1:50], vprofile(round(length(vprofile)/2), indices),'r');
legend('Womersley at R=0', 'UIV at R=0');
figfilename=['phantom_a_Delta' int2str(delta) '_WomersleyComparison_R=0b.png'];
print( '-painters', '-dpng', '-r600', figfilename);

womerpro=real(vel3);
UIVpro=interp1(r(fliplr(1:round(length(vprofile)/2))),vprofile_average(:,indices),[0:0.01:1].*R(indices(t))*1000);
meanWomerror=mean(mean(abs(womerpro-UIVpro')))
stdWomerror=std(mean(abs(womerpro-UIVpro')))
percenterror=meanWomerror./(max(max(womerpro))).*100
stdpercenterror=stdWomerror./(max(max(womerpro))).*100
percenterrorAxis=abs(womerpro(:,1)-UIVpro(1,:)')./abs(womerpro(:,1));

%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;
heightBox=0.2;
widthBox=3;
[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
    pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    for n=1:length(vmag_smaller)
        v=round(vmag_smaller(n)*100);
        pos=find(velocities==v);
        intensity(pos,T)=intensity(pos,T)+1;
    end
    averageUIV(t)=mean(vmag_smaller);
    T=T+1;
end
intensity=flipud(intensity);

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/a/Doppler/a_Doppler.mat', vars{:});
[v_long,t_long]=size(average);
velocitiesDS=([0:v_long]-unwrap).*ourVpix.*0.01;
timesDS=time;

%spectra
fig=figure;
imagesc(timesDS, velocitiesDS, average); colormap gray;
intensity2= ind2rgb(flipud(intensity),jet);
hold on; g=imagesc(times, velocities./100,intensity2(:,indices,:));
set(g,'AlphaData',flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
set(gca, 'position', [0.2 0.2 0.8 0.65],'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.5 1.5]);
xlim([0 0.96]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
title({'Doppler spectrum   ',['overlayed with UIV spectrum \Deltat=' int2str(delta) '\tau    ']},'FontSize',8,'FontWeight','normal');

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5.5 5]);
set(gcf, 'PaperPosition', [0 0 5.5 5]);
figfilename=['phantom_a_Delta' int2str(delta) '_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%error
background=mean(average(1,:));
sumAverageDS=(velocitiesDS(1:end-1)'*ones(1,length(timesDS))).*(average-background);
weightedAverageDS=sum(sumAverageDS)./sum(average-background);
weightedAverageDS2=interp1(timesDS,weightedAverageDS,times);
meanDiffaveDS=mean(abs(averageUIV(indices)-weightedAverageDS2))
percentDiffaveDS=meanDiffaveDS/max(averageUIV).*100;

peak=max(average);
for t=1:length(average)
    threshDS=find(average(:,t)>(peak(t)/2));
    maxDS(t)=velocitiesDS(threshDS(end));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices)-maxDS2))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100
stdDiffmaxDS=std(abs(averageUIV(indices)-maxDS2))
stdpercentDiffmaxDS=stdDiffmaxDS/max(averageUIV).*100
%% -----1 velocity profile---------
yfrom=1; yto=241;
idx1=21;

for t=indices(idx1)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(102);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx1))/iastep2);
Lower=round(mean(lowerEdge(:,idx1))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);

hold on;
if delta==7
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
elseif delta==15
    plot(r,vprofile(Upper:Lower),'m','LineWidth',0.5);
elseif delta==31
    plot(r,vprofile(Upper:Lower),'r','LineWidth',0.5);
elseif delta == 47
    plot(r,vprofile(Upper:Lower),'Color',[1 0.5 0],'LineWidth',0.5);
elseif delta==63
    plot(r,vprofile(Upper:Lower),'g','LineWidth',0.5);
elseif delta==79
    plot(r,vprofile(Upper:Lower),'c','LineWidth',0.5);
elseif delta==95
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
end

%% -----1 velocity profile---------
yfrom=1; yto=241;
idx2=28;

for t=indices(idx2)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(103);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx2))/iastep2);
Lower=round(mean(lowerEdge(:,idx2))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);

hold on;
if delta==7
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
elseif delta==15
    plot(r,vprofile(Upper:Lower),'m','LineWidth',0.5);
elseif delta==31
    plot(r,vprofile(Upper:Lower),'r','LineWidth',0.5);
elseif delta == 47
    plot(r,vprofile(Upper:Lower),'Color',[1 0.5 0],'LineWidth',0.5);
elseif delta==63
    plot(r,vprofile(Upper:Lower),'g','LineWidth',0.5);
elseif delta==79
    plot(r,vprofile(Upper:Lower),'c','LineWidth',0.5);
elseif delta==95
    plot(r,vprofile(Upper:Lower),'Color',[0.55 0 0.55],'LineWidth',0.5);
end
%% ------write out errors file
filename=['phantom_a_Delta' int2str(delta) 'errors.txt'];

fid=fopen(filename,'w');
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffmaxDS,percentDiffmaxDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = '%5.3f, %5.3f \n\n';
fprintf(fid,formatSpec,meanWomerror,percenterror);

formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,stdDiffmaxDS,stdpercentDiffmaxDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,stdFlowerror,stdpercentFlowerror);
formatSpec = '%5.3f, %5.3f \n\n';
fprintf(fid,formatSpec,stdWomerror,stdpercenterror);
fclose('all');
end
%%---conventional----
%% -----set up------
%clear VxC2 VyC2 VxC3 VyC3 VxF VyF
name=['../data/flow_phantom/a/UIV/2cm/conventional/phantom_a_128lines_conventional_20000depth7000000freq_out']
load([name '.mat']);

lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(lines, depth, Vx, Vy, 0,'glyce',period, split);

indices=[21:50 1:20];
centre=[1.9,1.1];

%--sort out the time steps--
timestep=(lines)*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% -----velocity profile---------
clear vprofile
yfrom=100; yto=180;
for t=indices
    vprofile(:,t)=mean(-VxC([1:22 36:63],yfrom:yto,t),1);
end
r=([1:length(vprofile(:,t))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;

%% ----flow rate-comparison----
%---flow rate UIV---
[flow2, Flow, upperEdge, lowerEdge]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 390, 440, 670, 730);
flow=flow2(indices);
figure(100); hold on;
plot(times,flow,'b','LineWidth',0.1);

set(gca,'Position',[0.2 0.2 0.76 0.76]);
Tindices=[30:length(flowT) 1:29];
hold on; plot(timeT, flowT(Tindices),'k','LineWidth',0.1);
hold on; plot(timeT,flowTMax(Tindices),'k:','LineWidth',0.1);
hold on; plot(timeT,flowTMin(Tindices),'k:','LineWidth',0.1);
xlabel('time /s','FontSize',8);
ylabel('flow rate /ml min^{-1}','FontSize',8);
strings={'\Delta t=15','\Delta t=31','\Delta t=47','\Delta t=63','\Delta t=79','\Delta t=95','conventional',['flow meter' 10 'mean+/-std']};
legend(strings,'Position',[0.72 0.65 0.2 0.2],'FontSize',6);
legend('boxoff')
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['phantom_a_flow_rate_comparison.png'];
print( '-painters', '-dpng', '-r600', figfilename);

flowT2 = interp1([timeT], [flowT(Tindices)],times(1:end-1));
meanFlowerror=mean(abs(flow(1:end-1)-flowT2'))
percentFlowerror=meanFlowerror./(max(flow)).*100

%% ----womersley profile----
figure;
for t=1:50
    R(t)=ypix*mean(lowerEdge(:,t)-upperEdge(:,t))/2;
end
Rmean=mean(R);
flow_f=flowT(Tindices);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));

vmean=flowT(Tindices)./(pi*Rmean^2);
vmean2=vmean./(60*1000*1000);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
vel3= womer_Vmean(timeT(pos), 20, vmean2(pos), Rmean,'old_e',nu);

vprofile_average=(vprofile(round(length(vprofile)/2):end,:)+vprofile(fliplr(1:round(length(vprofile)/2)),:))/2;

figure;  plot([1:50], real(vel3(:,1)),'k');
hold on; plot([1:50], vprofile(round(length(vprofile)/2), indices),'r');
legend('Womersley at R=0', 'UIV at R=0');
figfilename=['phantom_a_conventional_WomersleyComparison_R=0b.png'];
print( '-painters', '-dpng', '-r600', figfilename);

womerpro=real(vel3);
UIVpro=interp1(r(fliplr(1:round(length(vprofile)/2))),vprofile_average(:,indices),[0:0.01:1].*R(indices(t))*1000);
meanWomerror=mean(mean(abs(womerpro-UIVpro')))
percenterror=meanWomerror./(max(max(womerpro))).*100
percenterrorAxis=abs(womerpro(:,1)-UIVpro(1,:)')./abs(womerpro(:,1));

%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;
heightBox=0.2;
widthBox=3;
[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
    pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    for n=1:length(vmag_smaller)
        v=round(vmag_smaller(n)*100);
        pos=find(velocities==v);
        intensity(pos,T)=intensity(pos,T)+1;
    end
    averageUIV(t)=mean(vmag_smaller);
    T=T+1;
end
intensity=flipud(intensity);

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/a/Doppler/a_Doppler.mat', vars{:});
[v_long,t_long]=size(average);
velocitiesDS=([0:v_long]-unwrap).*ourVpix.*0.01;
timesDS=time;

%spectra
fig=figure;
imagesc(timesDS, velocitiesDS, average); colormap gray;
intensity2= ind2rgb(flipud(intensity),jet);
hold on; g=imagesc(times, velocities./100,intensity2(:,indices,:));
set(g,'AlphaData',flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
set(gca, 'position', [0.2 0.2 0.8 0.65],'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.5 1.5]);
xlim([0 0.96]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
title({'Doppler spectrum   ','overlayed with UIV spectrum conventional              '},'FontSize',8,'FontWeight','normal');

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5.5 5]);
set(gcf, 'PaperPosition', [0 0 5.5 5]);
figfilename=['phantom_a_conventional_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%error
background=mean(average(1,:));
sumAverageDS=(velocitiesDS(1:end-1)'*ones(1,length(timesDS))).*(average-background);
weightedAverageDS=sum(sumAverageDS)./sum(average-background);
weightedAverageDS2=interp1(timesDS,weightedAverageDS,times);
meanDiffaveDS=mean(abs(averageUIV(indices)-weightedAverageDS2))
percentDiffaveDS=meanDiffaveDS/max(averageUIV).*100;

peak=max(average);
for t=1:length(average)
    threshDS=find(average(:,t)>(peak(t)/2));
    maxDS(t)=velocitiesDS(threshDS(end));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices)-maxDS2))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100

%% -----1 velocity profile---------
yfrom=1; yto=241;
for t=indices(idx1)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(102);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx1))/iastep2);
Lower=round(mean(lowerEdge(:,idx1))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);

hold on;
plot(r,vprofile(Upper:Lower),'b','LineWidth',0.5);
xlim([0 3]);
ylim([0 1.2]);
set(gca,'FontSize',8);
title(['time = ' num2str(idx1/50.*times(end),2) ' s'],'FontSize',8);
xlabel('radial position /mm','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [4.1 5]);
set(gcf, 'PaperPosition', [0 0 4.1 5]);
figfilename=['phantom_a_vel_profile_comparison' int2str(idx1) '.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% -----1 velocity profile---------
yfrom=1; yto=241;
for t=indices(idx2)
    vmag=(VxC(:,:,t).^2+VyC(:,:,t).^2).^0.5;
    vmag=vmag.*sign(VxC(:,:,t));
    vprofile=mean(-vmag([11:45],yfrom:yto),1);
end
figure(103);
r=([1:length(vprofile)]).*iastep2.*ypix.*1;
Upper=round(mean(upperEdge(:,idx2))/iastep2);
Lower=round(mean(lowerEdge(:,idx2))/iastep2);
r=(r-ypix*iastep2*(Upper+Lower)/2)*1000;
r=r(Upper:Lower);

hold on;
plot(r,vprofile(Upper:Lower),'b','LineWidth',0.5);
xlim([0 3]);
ylim([0 1.2]);
set(gca,'FontSize',8);
title(['time = ' num2str(idx2/50.*times(end),2) ' s'],'FontSize',8);
xlabel('radial position /mm','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [4.1 5]);
set(gcf, 'PaperPosition', [0 0 4.1 5]);
figfilename=['phantom_a_vel_profile_comparison' int2str(idx2) '.png'];
print( '-painters', '-dpng', '-r600', figfilename);


%% ------write out errors file
filename=['phantom_a_Delta_normal_errors.txt'];
fid=fopen(filename,'w');
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffaveDS,percentDiffaveDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanDiffmaxDS,percentDiffmaxDS);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = '%5.3f, %5.3f \n';
fprintf(fid,formatSpec,meanWomerror,percenterror);
fclose('all');