From a296dcf2a6e2050911b77ca8d0f7acccacf44024 Mon Sep 17 00:00:00 2001 From: "Bhattacharya, Sayantan" Date: Fri, 30 Apr 2021 04:02:37 -0400 Subject: [PATCH] Add files via upload adding postprocessing codes for MC uncertainty --- postprocess_codes/Jet_data_postprocess.m | 839 +++++++++ .../Jet_data_postprocess_pranaonly.m | 569 ++++++ .../RMS_line_profile_for_2003B.m | 134 ++ postprocess_codes/RMS_spatial_profiles.m | 456 +++++ .../RMS_spatial_profiles_2Dmap.m | 420 +++++ postprocess_codes/compute_error_uncertainty.m | 783 +++++++++ .../convert_sim_images_to_davis.m | 35 + postprocess_codes/gradient_compact_rich.m | 99 ++ postprocess_codes/gradient_compactrich.m | 133 ++ .../gradient_dudxdvdy_compact_rich.m | 129 ++ postprocess_codes/in_out_directory_set.m | 79 + .../in_out_directory_set_03_09_2017.m | 80 + .../in_out_directory_set_03_09_2017_ws64.m | 79 + .../in_out_directory_set_07_30_2018_ws1.m | 79 + .../in_out_directory_set_07_30_2018_ws2.m | 80 + .../in_out_directory_set_10_22_2015.m | 79 + .../pivchal05b_ixx_postprocess.m | 681 ++++++++ postprocess_codes/plot_MC_uncertainty.m | 1534 +++++++++++++++++ .../plotting_jetdata_uncertainty.m | 51 + postprocess_codes/postprocess_prana_2003B.m | 602 +++++++ .../run_MC_uncertainty_analysis.m | 189 ++ postprocess_codes/socdiff.m | 56 + .../stagnation_flow_postprocess.m | 399 +++++ postprocess_codes/tridiagSolve.m | 31 + ...texring_original_uncertainty_postprocess.m | 599 +++++++ 25 files changed, 8215 insertions(+) create mode 100644 postprocess_codes/Jet_data_postprocess.m create mode 100644 postprocess_codes/Jet_data_postprocess_pranaonly.m create mode 100644 postprocess_codes/RMS_line_profile_for_2003B.m create mode 100644 postprocess_codes/RMS_spatial_profiles.m create mode 100644 postprocess_codes/RMS_spatial_profiles_2Dmap.m create mode 100644 postprocess_codes/compute_error_uncertainty.m create mode 100644 postprocess_codes/convert_sim_images_to_davis.m create mode 100644 postprocess_codes/gradient_compact_rich.m create mode 100644 postprocess_codes/gradient_compactrich.m create mode 100644 postprocess_codes/gradient_dudxdvdy_compact_rich.m create mode 100644 postprocess_codes/in_out_directory_set.m create mode 100644 postprocess_codes/in_out_directory_set_03_09_2017.m create mode 100644 postprocess_codes/in_out_directory_set_03_09_2017_ws64.m create mode 100644 postprocess_codes/in_out_directory_set_07_30_2018_ws1.m create mode 100644 postprocess_codes/in_out_directory_set_07_30_2018_ws2.m create mode 100644 postprocess_codes/in_out_directory_set_10_22_2015.m create mode 100644 postprocess_codes/pivchal05b_ixx_postprocess.m create mode 100644 postprocess_codes/plot_MC_uncertainty.m create mode 100644 postprocess_codes/plotting_jetdata_uncertainty.m create mode 100644 postprocess_codes/postprocess_prana_2003B.m create mode 100644 postprocess_codes/run_MC_uncertainty_analysis.m create mode 100644 postprocess_codes/socdiff.m create mode 100644 postprocess_codes/stagnation_flow_postprocess.m create mode 100644 postprocess_codes/tridiagSolve.m create mode 100644 postprocess_codes/vortexring_original_uncertainty_postprocess.m diff --git a/postprocess_codes/Jet_data_postprocess.m b/postprocess_codes/Jet_data_postprocess.m new file mode 100644 index 0000000..4a3084e --- /dev/null +++ b/postprocess_codes/Jet_data_postprocess.m @@ -0,0 +1,839 @@ +% Comapre Prana and Davis Solutions for experimental Jet data +clear all +strcase='cropped'; +printfig=0; +plotfig=1; +savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\'; +N=496; +truesoldir='Z:\Jet_Uncertainty_Data\HDR_data_F001\'; +% davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303_y_188_frame_03_501\TR_PIV_MP(3x16x16_75%ov)_01\'; % without vector pre and post processing +% davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303_y_188_frame_03_501\TR_PIV_MP(1x32x32_87%ov)(3x16x16_75%ov)_nopostprocess\'; % without vector pre and post processing +davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501\TR_PIV_MP(1x32x32_75%ov)(3x16x16_75%ov)_postprocess\'; % without vector pre and post processing + +% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\'; +% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\newixx\'; +% soldir2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\'; + +% soldir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\'; +% soldir2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\ws32IM\'; + +soldir='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Jetdata\'; +soldir2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Jetdata\'; + + +% soldir2='Z:\Planar_Uncertainty_work\Results\Experimental_Jet\unsteady_inviscid_core\New_processing_test_compare_prana_and_Davis\cropped_x_303_y_188\set3_500_firstpass64\';% Image Matching +soldircs='Z:\Planar_Uncertainty_work\Results\Experimental_Jet\unsteady_inviscid_core\New_processing_test_compare_prana_and_Davis\cropped_x_303_y_188\Correlation_Statistics\with_filter\'; % Corr stat + +% outputdir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\Jetdata\newixx2\'; +outputdir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\Jetdata\ws32IM\'; + + +biasfact=0; + +if ~exist(outputdir,'dir') + mkdir(outputdir); +end + +addpath Y:\Projects\TurbulentFlame\analysis\src\readimx_v2.0_win64\readimx; +basename1='PIV_scc_jet_cropped_pass4_'; +basename2='B'; +basename3='PIV_scc_jet_cropped_pass4_'; +% c1=4;c2=26;%303+14-1 303+102-1 +% r1=4; r2=36; %188+14-1 188+142-1 +c1=5;c2=23;%303+14-1 303+102-1 +r1=5; r2=35; %188+14-1 188+142-1 + +lw=3;fs=20; + +I=78;J=70; +% Sx=66;Sy=46; + +% Sx=40;Sy=30; +Sx=40;Sy=25; + +for i=1:N + + truesol=importdata([truesoldir,'B',num2str(i+2,'%05.0f'),'.dat']); + sccsol=load([soldir,basename1,num2str(i+2,'%05.0f'),'.mat']); + sccsolcs=load([soldircs,basename3,num2str(i+2,'%05.0f'),'.mat']); + sccsol2=load([soldir2,basename3,num2str(i+2,'%05.0f'),'.mat']); + %% Davis Result + M=readimx([davisdir,basename2,num2str(i,'%05.0f'),'.vc7']); + + if i==1 + [D] = create2DVec(M.Frames{1,1}); + Xd=D.X; + Yd=D.Y; + Xd=Xd'; + Yd=Yd'; + if strcmpi(strcase,'whole') + Xd=Xd(r1:r2,c1:c2); + Yd=Yd(r1:r2,c1:c2); + elseif strcmpi(strcase,'cropped') + Xd=Xd(r1:r2,c1:c2)+302; + Yd=Yd(r1:r2,c1:c2)+187; + end + + clear D; + end + + + Udtemp=double(cell2mat(M.Frames{1,1}.Components{1,1}.Planes)'); + Vdtemp=double(cell2mat(M.Frames{1,1}.Components{2,1}.Planes)'); + Ud(:,:,i)=Udtemp(r1:r2,c1:c2); + Vd(:,:,i)=Vdtemp(r1:r2,c1:c2); + %% Get Correlation Statistics estimates + + Unrxtemp=cell2mat(M.Frames{1,1}.Components{22,1}.Planes)'; + Unrytemp=cell2mat(M.Frames{1,1}.Components{23,1}.Planes)'; + Unbxtemp=cell2mat(M.Frames{1,1}.Components{25,1}.Planes)'; + Unbytemp=cell2mat(M.Frames{1,1}.Components{26,1}.Planes)'; + Untrtemp=cell2mat(M.Frames{1,1}.Components{21,1}.Planes)'; + Untbtemp=cell2mat(M.Frames{1,1}.Components{24,1}.Planes)'; + + Unrx(:,:,i)=Unrxtemp(r1:r2,c1:c2); + Unry(:,:,i)=Unrytemp(r1:r2,c1:c2); + Unbx(:,:,i)=Unbxtemp(r1:r2,c1:c2); + Unby(:,:,i)=Unbytemp(r1:r2,c1:c2); + Untr(:,:,i)=Untrtemp(r1:r2,c1:c2); + Untb(:,:,i)=Untbtemp(r1:r2,c1:c2); + %% True solution + if i==1 + X=truesol.data(:,1); + Y=truesol.data(:,2); + flag=truesol.data(:,6); + X=reshape(X,J,I); + Y=reshape(Y,J,I); + X=X';Y=Y'; + +% tsx=6:2:71; +% tsy=5:2:50; + tsx=8:2:69; + tsy=7:2:43; + + X=X(tsx,tsy); + Y=Y(tsx,tsy); + + Xh=(Xd-105)./(10.2*7.2); + Yh=-((Yd-254)./(10.2*7.2)); + end + + U=truesol.data(:,3); + V=truesol.data(:,4); + Ut=reshape(U,J,I)'; + Vt=reshape(V,J,I)'; + + + Utemp=Ut(tsx,tsy); + Vtemp=Vt(tsx,tsy); + + + + + Utrue(:,:,i)=Utemp; + Vtrue(:,:,i)=Vtemp; + + %% prana Solution + + if i==1 + + Xs=sccsol.X+302; + Ys=sccsol.Y+187; + + % + Xs=Xs(r1:r2,c1:c2); + Ys=Ys(r1:r2,c1:c2); + end + Ustemp=sccsol.U(end:-1:1,:); + Vstemp=-sccsol.V(end:-1:1,:); + Us(:,:,i)=Ustemp(r1:r2,c1:c2); + Vs(:,:,i)=Vstemp(r1:r2,c1:c2); + + + %% Get IXX/IYY estimates + + % Storing Ixx and other parameters. + tppr=reshape(sccsol.uncertainty(:,1),Sx,Sy); + tppr=tppr(end:-1:1,:); + PPR(:,:,i)=tppr(r1:r2,c1:c2); + + tmi=reshape(sccsol.uncertainty(:,5),Sx,Sy); + tmi=tmi(end:-1:1,:); + MI(:,:,i)=tmi(r1:r2,c1:c2); + + Autod=reshape(sccsol.uncertainty(:,6),Sx,Sy); + Autod=Autod(end:-1:1,:); + Autod=Autod(r1:r2,c1:c2); + + Ixx=reshape(sccsol.uncertainty(:,7),Sx,Sy); + Ixx=Ixx(end:-1:1,:); + Ixx=Ixx(r1:r2,c1:c2); + + Iyy=reshape(sccsol.uncertainty(:,8),Sx,Sy); + Iyy=Iyy(end:-1:1,:); + Iyy=Iyy(r1:r2,c1:c2); + + %Gcc peak location + Gpx=reshape(sccsol.uncertainty(:,15),Sx,Sy); + Gpx=Gpx(end:-1:1,:); + Gpx=Gpx(r1:r2,c1:c2); + + Gpy=reshape(sccsol.uncertainty(:,16),Sx,Sy); + Gpy=Gpy(end:-1:1,:); + Gpy=Gpy(r1:r2,c1:c2); + + %Scc peak location + Spx=reshape(sccsol.uncertainty(:,17),Sx,Sy); + Spx=Spx(end:-1:1,:); + Spx=biasfact*Spx(r1:r2,c1:c2); + + Spy=reshape(sccsol.uncertainty(:,18),Sx,Sy); + Spy=Spy(end:-1:1,:); + Spy=biasfact*Spy(r1:r2,c1:c2); + + %10:89,22:101; + %Bias error + mewx=abs(Spx-Gpx); + mewy=abs(Spy-Gpy); + % correct for gradient + [Udiff]=socdiff(Us(:,:,i),2,1); + [Vdiff]=socdiff(Vs(:,:,i),2,2); + + Ixx= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 )); + Iyy= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 )); + + %Scaling for ixx + Ixx=Ixx./sqrt(MI(:,:,i)); + Iyy=Iyy./sqrt(MI(:,:,i)); + +% Ixxtemp(:,:,i)=Ixx; +% Iyytemp(:,:,i)=Iyy; + + % Final second order moment uncertainty + IXX1(:,:,i)=sqrt(Ixx.^2+mewx.^2); + IYY1(:,:,i)=sqrt(Iyy.^2+mewy.^2); + +% IXX1(:,:,i)=Ixx; +% IYY1(:,:,i)=Iyy; + + %% Get IMx/IMy estimates + % Image Matching uncertainty + timx=reshape(sccsol2.uncertainty(:,9),Sx,Sy); + timx=timx(end:-1:1,:); + Imx(:,:,i)=timx(r1:r2,c1:c2); + + timy=reshape(sccsol2.uncertainty(:,10),Sx,Sy); + timy=timy(end:-1:1,:); + Imy(:,:,i)=timy(r1:r2,c1:c2); + + tnim=reshape(sccsol2.uncertainty(:,11),Sx,Sy); + tnim=tnim(end:-1:1,:); + Nim(:,:,i)=tnim(r1:r2,c1:c2); + + %% Get prana correlation statistics estimates +% tcsx=reshape(sccsolcs.uncertainty(:,25),Sx,Sy); +% tcsx=tcsx(end:-1:1,:); +% Uncsx(:,:,i)=tcsx(r1:r2,c1:c2); +% +% tcsy=reshape(sccsolcs.uncertainty(:,26),Sx,Sy); +% tcsy=tcsy(end:-1:1,:); +% Uncsy(:,:,i)=tcsy(r1:r2,c1:c2); + + %% Error + error_Us(:,:,i)=(Utrue(:,:,i)-Us(:,:,i)); + error_Vs(:,:,i)=(Vtrue(:,:,i)-Vs(:,:,i)); + error_Ud(:,:,i)=(Utrue(:,:,i)-Ud(:,:,i)); + error_Vd(:,:,i)=(Vtrue(:,:,i)-Vd(:,:,i)); + UdUs(:,:,i)=(Ud(:,:,i)-Us(:,:,i)); + VdVs(:,:,i)=(Vd(:,:,i)-Vs(:,:,i)); + + %% Plot Velcity distributions +% % if printfig==1 +% % if i==1%>=1 || i<=10 +% % figure;pcolor(Xh,Yh,Utrue(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Utrue'); +% % print(gcf,'-dpng',[outputdir,'Utrue_frame_',num2str(i+2,'%05.0f'),'.png']); +% % figure;pcolor(Xh,Yh,Us(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Uprana'); +% % print(gcf,'-dpng',[outputdir,'Uprana_frame_',num2str(i+2,'%05.0f'),'.png']); +% % figure;pcolor(Xh,Yh,Ud(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Udavis'); +% % print(gcf,'-dpng',[outputdir,'Udavis_frame_',num2str(i+2,'%05.0f'),'.png']); +% % +% % figure;pcolor(Xh,Yh,Vtrue(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vtrue'); +% % print(gcf,'-dpng',[outputdir,'Vtrue_frame_',num2str(i+2,'%05.0f'),'.png']); +% % figure;pcolor(Xh,Yh,Vs(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vprana'); +% % print(gcf,'-dpng',[outputdir,'Vprana_frame_',num2str(i+2,'%05.0f'),'.png']); +% % figure;pcolor(Xh,Yh,Vd(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vdavis'); +% % print(gcf,'-dpng',[outputdir,'Vdavis_frame_',num2str(i+2,'%05.0f'),'.png']); +% % +% % +% % figure;%set(gcf,'DefaultLineLineWidth',1); +% % set(gca,'FontSize',10);hold on; +% % quiver(Xh,Yh,Utrue(:,:,i),Vtrue(:,:,i),1,'color','red');axis image; +% % quiver(Xh,Yh,Us(:,:,i),Vs(:,:,i),1,'color','black'); +% % quiver(Xh,Yh,Ud(:,:,i),Vd(:,:,i),1,'color','blue'); +% % hold off; +% % box on; +% % xlabel('x/h');ylabel('y/h'); +% % title(['Velocity field at frame',num2str(i+2,'%05.0f')]); +% % % legend('True','prana','davis'); +% % print(gcf,'-dpng',[outputdir,'Vectorfield_frame2_',num2str(i+2,'%05.0f'),'.png'],'-r1200'); +% % +% % % keyboard; +% % +% % pause(0.1); +% % close all; +% % pause(0.2); +% % end +% % end + % if i==1 +% figure; +% subplot(1,2,2);pcolor(Xh,Yh,error_V(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Vtrue');colorbar; +% subplot(1,2,1);pcolor(Xh,Yh,error_U(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Utrue');colorbar; +% if printfig==1 +% set(gcf,'units','normalized','outerposition',[0 0 1 1]); +% print(gcf,'-dpng',[outputdir,'Error_frame_',num2str(i+2,'%05.0f'),'.png'],'-r300'); +% end +% end + +% keyboard; + if rem(i,100)==0 + i + end +end +%% Save error and uncertainties +% save([savematdir,'Jet_uncertainty_data_x_303_y_188_firstpass64_cclsq.mat'],'error_Us','error_Vs','error_Ud','error_Vd','IXX1','IYY1','Imx','Imy','Unrx','Unry','Unbx','Unby'); + +% keyboard; +%% +teUs=error_Us(:); +% teUs(abs(teUs)>0.5)=[]; +teUd=error_Ud(:); +% teUd(abs(teUd)>0.5)=[]; +teUsd=UdUs(:); +% teUsd(abs(teUsd)>0.5)=[]; + + +teVs=error_Vs(:); +% teVs(abs(teVs)>0.5)=[]; +teVd=error_Vd(:); +% teVd(abs(teVd)>0.5)=[]; +teVsd=VdVs(:); +% teVsd(abs(teVsd)>0.5)=[]; + +vec=linspace(-0.5,0.5,50); + +Neus=histc(teUs,vec); +Neud=histc(teUd,vec); +Neusd=histc(teUsd,vec); + +Nevs=histc(teVs,vec); +Nevd=histc(teVd,vec); +Nevsd=histc(teVsd,vec); + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); +hold on; +plot(vec,Neus,'b-'); +plot(vec,Neud,'r-'); +% plot(vec,Neusd,'k-'); +% legend('Utrue-Uprana','Utrue-Udavis','Udavis-Uprana'); +legend('Utrue-Uprana','Utrue-Udavis'); +title('U error histogram'); +hold off; +grid on; +% keyboard; +if printfig==1 +print(gcf,'-dpng',[outputdir,'UerrorHist2','.png']); +end +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); +hold on; +plot(vec,Nevs,'b-'); +plot(vec,Nevd,'r-'); +% plot(vec,Nevsd,'k-'); +% legend('Vtrue-Vprana','Vtrue-Vdavis','Vdavis-Vprana'); +legend('Vtrue-Vprana','Vtrue-Vdavis'); +title('V error histogram'); +hold off; +grid on; +if printfig==1 + print(gcf,'-dpng',[outputdir,'VerrorHist2','.png']); +end +% keyboard; +%% Instantaneous U Velocity profile + xloc=16;yloc=14; %(378,254, 376,253) +% xloc=21;yloc=14; %(378,254, 376,253) + tstamp=10; + Uinsttrue=squeeze(Utrue(:,xloc,tstamp)); + Uinstsol=squeeze(Us(:,xloc,tstamp)); + Uinstdav=squeeze(Ud(:,xloc,tstamp)); + figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + plot(Yh(:,1),Uinsttrue,'k-');hold on; + plot(Yh(:,1),Uinstsol,'b-'); + plot(Yh(:,1),Uinstdav,'r-');hold off; + axis([-1.4 1.4 -0.4 4.2]); + xlabel('Yh'); + ylabel('U'); + title(['Instantaneous U Velocity at t=',num2str(0.1*tstamp,'%02.1f'),' ms']); + legend('True','prana','davis'); + box on; +% set(gcf,'units','normalized','outerposition',[0 0 1 1]); +if printfig==1 + print(gcf,'-dpng',[outputdir,'InstantU.png']); +end + %% velocity time series + + UHDts=squeeze(Utrue(yloc,xloc,:)); + UMSts=squeeze(Us(yloc,xloc,:)); + UDVts=squeeze(Ud(yloc,xloc,:)); + + ts=0:0.1:0.1*length(UHDts)-0.1; + ts=ts'; + figure;set(gca,'FontSize',fs); + hold on; + plot(ts,UHDts,'k-','LineWidth',3); + plot(ts,UMSts,'b-','LineWidth',1.5); + plot(ts,UDVts,'r-','LineWidth',1.5); + axis([0 40 2.95 4.05]); + hold off; + xlabel('t (ms)'); + ylabel('U (px)'); + title(['Velocity time series at X=',num2str(X(1,xloc),'%02.1f'),'Y=',num2str(Y(yloc,1),'%02.1f')]); + legend('True','Prana','Davis'); + set(gcf,'units','normalized','outerposition',[0 0 1 1]); + box on; + if printfig==1 + print(gcf,'-dpng',[outputdir,'Velocity time series.png']); + end + %% Error Time Series + eUts=squeeze(error_Us(yloc,xloc,:)); + eUdts=squeeze(error_Ud(yloc,xloc,:)); + eVts=squeeze(error_Vs(yloc,xloc,:)); + eMts=sqrt(eUts.^2+eVts.^2); + meUts=mean(eUts).*ones(length(eUts)); + meUdts=mean(eUdts).*ones(length(eUdts)); + + figure;set(gca,'FontSize',fs); + plot(ts,eUts,'b-');hold on; + plot(ts,meUts,'b--','Linewidth',2);hold off; + axis([0 max(ts(:)) -0.19 0.19]); + xlabel('t (ms)'); + ylabel('e_U prana(px)'); + box on; + st1='Error U time series (prana)'; + st4='Error U time series (Davis)'; + st2=['At X=',num2str(X(1,xloc),'%02.1f'),',Y=',num2str(Y(yloc,1),'%02.1f')]; + st3=['mean error=',num2str(meUts(1),'%02.3f')]; + st5=['mean error=',num2str(meUdts(1),'%02.3f')]; + + title([st1,'\newline',st2,'\newline',st3]); + if printfig==1 + print(gcf,'-dpng',[outputdir,'Error U (Prana) time series.png']); + end + + + figure;set(gca,'FontSize',fs); + plot(ts,eUdts,'r-');hold on; + plot(ts,meUdts,'r--','Linewidth',2);hold off; + axis([0 max(ts(:)) -0.19 0.19]);box on; + xlabel('t (ms)'); + ylabel('e_U davis(px)'); + title([st4,'\newline',st2,'\newline',st5]); + if printfig==1 + print(gcf,'-dpng',[outputdir,'Error U (Davis) time series.png']); + end + + %% Instantaneous Velocity and Error profile + + etemp1=squeeze(error_Us(:,xloc,:)); + rmserrU=rms(etemp1,2); + etemp5=squeeze(error_Ud(:,xloc,:)); + rmserrUd=rms(etemp5,2); + + etemp2=squeeze(IXX1(:,xloc,:)); + rmserrIxx=rms(etemp2,2); + etemp3=squeeze(Imx(:,xloc,:)); + rmserrImx=rms(etemp3,2); + etemp6=squeeze(Unrx(:,xloc,:)); + rmserrIxxd1=rms(etemp6,2); + etemp7=squeeze(Unbx(:,xloc,:)); + rmserrIxxd2=rms(etemp7,2); + rmserrIxxd=sqrt(rmserrIxxd1.^2+rmserrIxxd2.^2); + +% etemp8=squeeze(Uncsx(:,xloc,:)); +% rmserrUncsx=rms(etemp8,2); + + figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(Yh(:,1),rmserrU,'r-'); + plot(Yh(:,1),rmserrUd,'r--'); + plot(Yh(:,1),rmserrIxx,'k-'); + plot(Yh(:,1),rmserrImx,'c-'); + plot(Yh(:,1),rmserrIxxd,'m--'); +% plot(Yh(:,1),rmserrUncsx,'m-'); + hold off;box on; + axis([-1.4 1.4 0 0.2]); + xlabel('Yh'); + ylabel('RMS error, uncertainty'); + title(['RMS error over time at X=',num2str(X(1,xloc),'%02.1f')]); +% legend('rms e_u (Prana)','rms e_u (Davis)','rms IXX (Prana)','rms PDx(Prana)','rms UCSx(Davis)','rms UCSx (Prana)'); + set(gcf,'units','normalized','outerposition',[0 0 1 1]); + if printfig==1 + print(gcf,'-dpng',[outputdir,'RMS error2.png']); + end + +keyboard; + + + + + %% + if plotfig==1 + + error_U=abs(error_Us); + error_V=abs(error_Vs); + %% convert to vectors + IXX1=IXX1(:); + IYY1=IYY1(:); + Imx=Imx(:); + Imy=Imy(:); + error_U=error_U(:); + error_V=error_V(:); + + IXX2=Imx; + IYY2=Imy; + % keyboard; + % save('VRSOMnew_toterr.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2'); + % keyboard; + %% zero value removal + inx=(find(IXX1(:)==0)); + iny=(find(IYY1(:)==0)); + %keyboard; + IXX1(inx)=[]; + IYY1(iny)=[]; + Imx(inx)=[]; + Imy(iny)=[]; + error_U(inx)=[]; + error_V(iny)=[]; + %% + %nan value removal + inx2=(find(isnan(Imx))); + iny2=(find(isnan(Imy))); + + %keyboard; + IXX1(inx2)=[]; + IYY1(iny2)=[]; + Imx(inx2)=[]; + Imy(iny2)=[]; + error_U(inx2)=[]; + error_V(iny2)=[]; + + inx3=(find(isnan(IXX1))); + iny3=(find(isnan(IXX1))); + IXX1(inx3)=[]; + IYY1(iny3)=[]; + Imx(inx3)=[]; + Imy(iny3)=[]; + error_U(inx3)=[]; + error_V(iny3)=[]; + + [length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)] + + %% invalid vector removal + inx1=(find(error_U(:)>1)); + iny1=(find(error_V(:)>1)); +% inx1=(find(error_U(:)>0.5)); +% iny1=(find(error_V(:)>0.5)); +%keyboard; + IXX1(inx1)=[]; + IYY1(iny1)=[]; + Imx(inx1)=[]; + Imy(iny1)=[]; + error_U(inx1)=[]; + error_V(iny1)=[]; + + %% Find Coverage + + gmx=(find(error_U(:)<=1)); + gmy=(find(error_V(:)<=1)); + ngmx=length(gmx); + ngmy=length(gmy); + + %keyboard; + % N=r*s*t; + % [(N-ngmx)*100/N (N-ngmy)*100/N]; + + cnt=0; + for j=1:length(IXX1) + if error_U(j)<=IXX1(j) + cnt=cnt+1; + end + end + covx=100*cnt/ngmx; + + cnt=0; + for j=1:length(IYY1) + if error_V(j)<=IYY1(j) + cnt=cnt+1; + end + end + covy=100*cnt/ngmy; + + cnt=0; + for j=1:length(Imx) + if error_U(j)<=Imx(j) + cnt=cnt+1; + end + end + covx2=100*cnt/ngmx; + + cnt=0; + for j=1:length(Imy) + if error_V(j)<=Imy(j) + cnt=cnt+1; + end + end + covy2=100*cnt/ngmy; + + [covx covx2] + [covy covy2] + + % [rms(error_U) rms(error_V)] + + % IXX2=Imx; + % IYY2=Imy; + + % save('VRSOM.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2'); + + %% Plot error and uncertainty histogram + thresh=0.5; + + vec=linspace(0,thresh,60); + % vec2=linspace(0,0.05,100); + % vec3=linspace(0.03,0.1,100); + % vec4=linspace(0,0.03,100); + % vec5=linspace(0.06,0.1,100); + + Nrex=histc(error_U,vec); + Nrey=histc(error_V,vec); + Nrdx=histc(IXX1,vec); + Nrdy=histc(IYY1,vec); + Nrdx2=histc(Imx,vec); + Nrdy2=histc(Imy,vec); + + % Npprxl=histc(upprx(:,2),vec2); + % Npprxh=histc(upprx(:,1),vec3); + % + % Nppryl=histc(uppry(:,2),vec2); + % Nppryh=histc(uppry(:,1),vec3); + % + % Nmixl=histc(umix(:,2),vec4); + % Nmixh=histc(umix(:,1),vec5); + % + % Nmiyl=histc(umiy(:,2),vec4); + % Nmiyh=histc(umiy(:,1),vec5); + + ll=1:100:max(Nrex);%max(Nrdx); + ll2=1:100:max(Nrey);%max(Nrdy); + + + sigu=rms(error_U(error_U<=thresh)).*ones(length(ll),1); + sigv=rms(error_V(error_V<=thresh)).*ones(length(ll2),1); + + % mixx1=median(IXX1(IXX1<=thresh)).*ones(length(ll),1); + % miyy1=median(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + % + % mixx2=median(Imx(Imx<=thresh)).*ones(length(ll),1); + % miyy2=median(Imy(Imy<=thresh)).*ones(length(ll2),1); + + mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1); + miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + + mixx2=rms(Imx(Imx<=thresh)).*ones(length(ll),1); + miyy2=rms(Imy(Imy<=thresh)).*ones(length(ll2),1); + + figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; + plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--'); + % plot(vec2,Npprxl,'color',rgb('DarkRed')); + % plot(vec3,Npprxh,'color',rgb('DarkRed')); + % plot(vec4,Nmixl,'color',rgb('orange')); + % plot(vec5,Nmixh,'color',rgb('orange')); + hold off; + %legend('|errorx|','IXX AMM','IXX DMM','rms|errorx|'); + %legend('|e_x|','I_X_X','\delta_x','\sigma_{|e_x|}','\sim{I_X_X}','\sim{\delta_x}'); + title('Error and Uncertainty histogram(X)'); + %title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)])); + ylabel('No. of Data Points');xlabel('pixels'); + axis([0 thresh 0 inf]); + if printfig==1 + print(gcf,'-dpng',[outputdir,'Uhist2.png']); + end + figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; + plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--'); + % plot(vec2,Nppryl,'color',rgb('DarkCyan')); + % plot(vec3,Nppryh,'color',rgb('DarkCyan')); + % plot(vec4,Nmiyl,'color',rgb('SkyBlue')); + % plot(vec5,Nmiyh,'color',rgb('SkyBlue')); + hold off; + %legend('|Errory|','IYY AMM','IYY DMM','rms|errory|'); + %legend('|e_y|','I_Y_Y','\delta_y','\sigma_{|e_y|}','\sim{I_Y_Y}','\sim{\delta_y}'); + title('Error and Uncertainty histogram(Y)'); + %title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)])); + ylabel('No. of Data Points');xlabel('pixels'); + axis([0 thresh 0 inf]); + if printfig==1 + print(gcf,'-dpng',[outputdir,'Vhist2.png']); + end + + [sigu(1) mixx1(1) mixx2(1)] + [sigv(1) miyy1(1) miyy2(1)] + %% + plotflag=1; + + if plotflag==1; + + IXX2=Imx; + IYY2=Imy; + err_sccu=error_U; + err_sccv=error_V; + + Nbin=15; + + maxixx=thresh; + minixx=0; + vu=linspace(minixx,maxixx,Nbin); + maxiyy=thresh; + miniyy=0; + vv=linspace(miniyy,maxiyy,Nbin); + + rmserru1=zeros(1,length(vu)); + rmserrv1=zeros(1,length(vv)); + getixx1=zeros(1,length(vu)); + getiyy1=zeros(1,length(vv)); + numcx1=zeros(1,length(vu)); + numcy1=zeros(1,length(vv)); + covxbin1=zeros(1,length(vu)); + covybin1=zeros(1,length(vv)); + + rmserru2=zeros(1,length(vu)); + rmserrv2=zeros(1,length(vv)); + getixx2=zeros(1,length(vu)); + getiyy2=zeros(1,length(vv)); + numcx2=zeros(1,length(vu)); + numcy2=zeros(1,length(vv)); + covxbin2=zeros(1,length(vu)); + covybin2=zeros(1,length(vv)); + + %p=1; + + for p=1:length(vu) + if p=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1=1 || i<=10 + figure;pcolor(Xh,Yh,Utrue(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Utrue'); +% print(gcf,'-dpng',[outputdir,'Utrue_frame_',num2str(i+2,'%05.0f'),'.png']); + figure;pcolor(Xh,Yh,Us(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Uprana'); +% print(gcf,'-dpng',[outputdir,'Uprana_frame_',num2str(i+2,'%05.0f'),'.png']); + + figure;pcolor(Xh,Yh,Vtrue(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vtrue'); +% print(gcf,'-dpng',[outputdir,'Vtrue_frame_',num2str(i+2,'%05.0f'),'.png']); + figure;pcolor(Xh,Yh,Vs(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vprana'); +% print(gcf,'-dpng',[outputdir,'Vprana_frame_',num2str(i+2,'%05.0f'),'.png']); + + figure;%set(gcf,'DefaultLineLineWidth',1); + set(gca,'FontSize',10);hold on; + quiver(Xh,Yh,Utrue(:,:,i),Vtrue(:,:,i),1,'color','red');axis image; + quiver(Xh,Yh,Us(:,:,i),Vs(:,:,i),1,'color','black'); + hold off; + box on; + xlabel('x/h');ylabel('y/h'); + title(['Velocity field at frame',num2str(i+2,'%05.0f')]); +% legend('True','prana','davis'); +% print(gcf,'-dpng',[outputdir,'Vectorfield_frame2_',num2str(i+2,'%05.0f'),'.png'],'-r1200'); + +% keyboard; + +% pause(0.1); +% close all; +% pause(0.2); +% end +% end +% if i==1 + figure; + subplot(1,2,2);pcolor(Xh,Yh,error_Vs(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Vtrue');colorbar; + subplot(1,2,1);pcolor(Xh,Yh,error_Us(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Utrue');colorbar; +% if printfig==1 +% set(gcf,'units','normalized','outerposition',[0 0 1 1]); +% print(gcf,'-dpng',[outputdir,'Error_frame_',num2str(i+2,'%05.0f'),'.png'],'-r300'); +% end +% end + + keyboard; + %} + keyboard; + if rem(i,100)==0 + i + end +end +%% Save error and uncertainties +% save([savematdir,'Jet_uncertainty_data_x_303_y_188_firstpass64_cclsq.mat'],'error_Us','error_Vs','error_Ud','error_Vd','IXX1','IYY1','Imx','Imy','Unrx','Unry','Unbx','Unby'); + +error_U=abs(error_Us); +error_V=abs(error_Vs); + +Ixxtemp=IXX1; +Iyytemp=IYY1; + + figure(1); + subplot(2,3,1);imagesc(mean(error_U,3));colorbar;caxis([0 0.2]); + subplot(2,3,2);imagesc(mean(IXX1,3));colorbar;caxis([0 0.2]); + subplot(2,3,3);imagesc(nanmean(Imx,3));colorbar;caxis([0 0.1]); + subplot(2,3,4);imagesc(mean(error_V,3));colorbar;caxis([0 0.1]); + subplot(2,3,5);imagesc(mean(IYY1,3));colorbar;caxis([0 0.1]); + subplot(2,3,6);imagesc(nanmean(Imy,3));colorbar;caxis([0 0.1]); +% pause(0.2); + +if printfig==1 + print(gcf,'-dpng',[outputdir,'Mean_profiles.png']); +end + +% keyboard; + + + + +%% convert to vectors +IXX1=IXX1(:); +IYY1=IYY1(:); +Imx=Imx(:); +Imy=Imy(:); +error_U=error_U(:); +error_V=error_V(:); +% IXX2=Imx; +% IYY2=Imy; +% keyboard; +% save('VRSOMnew_toterr.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2'); +% keyboard; +%% zero value removal +inx=(find(IXX1(:)==0)); +iny=(find(IYY1(:)==0)); +%keyboard; +IXX1(inx)=[]; +IYY1(iny)=[]; +Imx(inx)=[]; +Imy(iny)=[]; +error_U(inx)=[]; +error_V(iny)=[]; +%% +%nan value removal +inx2=(find(isnan(Imx))); +iny2=(find(isnan(Imy))); + +%keyboard; +IXX1(inx2)=[]; +IYY1(iny2)=[]; +Imx(inx2)=[]; +Imy(iny2)=[]; +error_U(inx2)=[]; +error_V(iny2)=[]; + +inx3=(find(isnan(IXX1))); +iny3=(find(isnan(IXX1))); +IXX1(inx3)=[]; +IYY1(iny3)=[]; +Imx(inx3)=[]; +Imy(iny3)=[]; +error_U(inx3)=[]; +error_V(iny3)=[]; + +[length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)] + +%% invalid vector removal +inx1=(find(error_U(:)>1)); +iny1=(find(error_V(:)>1)); +% inx1=(find(error_U(:)>0.4)); +% iny1=(find(error_V(:)>0.4)); +%keyboard; +IXX1(inx1)=[]; +IYY1(iny1)=[]; +Imx(inx1)=[]; +Imy(iny1)=[]; +error_U(inx1)=[]; +error_V(iny1)=[]; + +%% Find Coverage + +% gmx=(find(error_U(:)<=1)); +% gmy=(find(error_V(:)<=1)); +nx=length(error_U); +ny=length(error_V); + +%keyboard; +% N=r*s*t; +% [(N-ngmx)*100/N (N-ngmy)*100/N]; + +cnt=0; +for j=1:nx + if error_U(j)<=IXX1(j) + cnt=cnt+1; + end +end +covx=100*cnt/nx; + +cnt=0; +for j=1:ny + if error_V(j)<=IYY1(j) + cnt=cnt+1; + end +end +covy=100*cnt/ny; + +cnt=0; +for j=1:nx + if error_U(j)<=Imx(j) + cnt=cnt+1; + end +end +covx2=100*cnt/nx; + +cnt=0; +for j=1:ny + if error_V(j)<=Imy(j) + cnt=cnt+1; + end +end +covy2=100*cnt/ny; + +[covx covx2] +[covy covy2] + +% [rms(error_U) rms(error_V)] + +% IXX2=Imx; +% IYY2=Imy; + +% save('VRSOM.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2'); + +%% Plot error and uncertainty histogram +thresh=0.4; + +vec=linspace(0,thresh,60); + + +Nrex=histc(error_U,vec); +Nrey=histc(error_V,vec); +Nrdx=histc(IXX1,vec); +Nrdy=histc(IYY1,vec); +Nrdx2=histc(Imx,vec); +Nrdy2=histc(Imy,vec); + + +ll=1:100:max(Nrex);%max(Nrdx); +ll2=1:100:max(Nrey);%max(Nrdy); + + +sigu=rms(error_U(error_U<=thresh)).*ones(length(ll),1); +sigv=rms(error_V(error_V<=thresh)).*ones(length(ll2),1); + +mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1); +miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + +mixx2=rms(Imx(Imx<=thresh)).*ones(length(ll),1); +miyy2=rms(Imy(Imy<=thresh)).*ones(length(ll2),1); + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--'); +hold off; +title(sprintf(['Ixx=',num2str(covx),' IMx=',num2str(covx2)])); +%title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 thresh 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outputdir,'Uhist2.png']); +end +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--'); +hold off; +title('Error and Uncertainty histogram(Y)'); +title(sprintf(['Iyy=',num2str(covy),' IMy=',num2str(covy2)])); +%title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 thresh 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outputdir,'Vhist2.png']); +end +[sigu(1) mixx1(1) mixx2(1)] +[sigv(1) miyy1(1) miyy2(1)] + + +IXX2=Imx; +IYY2=Imy; +err_sccu=error_U; +err_sccv=error_V; + +Nbin=10; + +maxixx=0.2; +minixx=0; +vu=linspace(minixx,maxixx,Nbin); +maxiyy=0.2; +miniyy=0; +vv=linspace(miniyy,maxiyy,Nbin); + +rmserru1=zeros(1,length(vu)); +rmserrv1=zeros(1,length(vv)); +getixx1=zeros(1,length(vu)); +getiyy1=zeros(1,length(vv)); +numcx1=zeros(1,length(vu)); +numcy1=zeros(1,length(vv)); +covxbin1=zeros(1,length(vu)); +covybin1=zeros(1,length(vv)); + +rmserru2=zeros(1,length(vu)); +rmserrv2=zeros(1,length(vv)); +getixx2=zeros(1,length(vu)); +getiyy2=zeros(1,length(vv)); +numcx2=zeros(1,length(vu)); +numcy2=zeros(1,length(vv)); +covxbin2=zeros(1,length(vu)); +covybin2=zeros(1,length(vv)); + +%p=1; + +for p=1:length(vu) + if p=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1veccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; +%} + +% % %% Extract RMS profile at x=96th column +% % % RMS error +% % c1=96; +% % Upt=rms(squeeze(A1.err_up(end:-1:1,c1,:)),2); +% % Vpt=rms(squeeze(A1.err_vp(end:-1:1,c1,:)),2); +% % Velpt=sqrt(Upt.^2+Vpt.^2); +% % c2=95; +% % Udt=rms(squeeze(A1.err_ud(end:-1:1,c2,:)),2); +% % Vdt=rms(squeeze(A1.err_vd(end:-1:1,c2,:)),2); +% % Veldt=sqrt(Udt.^2+Vdt.^2); +% % +% % %RMS uncertainty +% % UMCt=rms(squeeze(A1.err_up(end:-1:1,c1,:)),2); +% % VMCt=rms(squeeze(A1.err_vp(end:-1:1,c1,:)),2); +% % VelMC=sqrt(Upt.^2+Vpt.^2); +% % +% % %% Plot RMS profiles +% % +% % figure;hold on; +% % plot(Velpt,1:63,'k-'); +% % plot(Veldt,1:62,'k--'); + +%% Extract RMS profile at x=96th column +% RMS error + + %% Eliminate bad vectors + veccutoff=0.1; +c1=96; +Upt=squeeze(A1.err_up(end:-1:1,c1,:));Upt(Upt>veccutoff)=nan; +Vpt=squeeze(A1.err_vp(end:-1:1,c1,:));Vpt(Vpt>veccutoff)=nan; +Velpt=sqrt(Upt.^2+Vpt.^2); +[inpy,inpx]=find(isnan(Velpt)); +Velptrms=sqrt(nanmean(Velpt.^2,2)); +c2=96; +Udt=squeeze(A1.err_ud(end:-1:1,c2,:));Udt(Udt>veccutoff)=nan; +Vdt=squeeze(A1.err_vd(end:-1:1,c2,:));Vdt(Vdt>veccutoff)=nan; +Veldt=sqrt(Udt.^2+Vdt.^2); +Veldtrms=sqrt(nanmean(Veldt.^2,2)); + +%RMS uncertainty +UMCt=rms(squeeze(A1.err_up(end:-1:1,c1,:)),2); +VMCt=rms(squeeze(A1.err_vp(end:-1:1,c1,:)),2); +VelMC=sqrt(Upt.^2+Vpt.^2); + + +%% Plot RMS profiles + +figure;hold on; +plot(Velptrms,1:63,'k-'); +plot(Veldtrms,1:62,'k--'); + diff --git a/postprocess_codes/RMS_spatial_profiles.m b/postprocess_codes/RMS_spatial_profiles.m new file mode 100644 index 0000000..a46bc09 --- /dev/null +++ b/postprocess_codes/RMS_spatial_profiles.m @@ -0,0 +1,456 @@ +%% Plot RMS spatial profiles for error and uncertainty +clear; +basedir='Z:\Planar_Uncertainty_work\Results\final_plots\07_31_2018\'; +% casename='WS64'; +% casename='WS32'; +% casename='WS32\New_gradient\'; +% casename='WS32\SOC_gradient\'; +casename='WS64\SOC_gradient\'; +outdir=[basedir,'spatial_RMS_profiles\',casename,filesep]; +if ~exist(outdir,'dir'); + mkdir(outdir); +end +savemat=1; +saveplot=1; +% keyboard; +%% 2003B +A1=load([basedir,casename,'\matfiles\PivChal03B.mat']); +% Upmean=mean(A1.Up,3); +% Vpmean=mean(A1.Vp,3); +% figure;imagesc(A1.Xp(:),A1.Yp(:),Upmean); +Xp=A1.Xp/max(A1.Xp(:)); +Yp=A1.Yp/max(A1.Yp(:)); +velmagmean=mean((A1.Up.^2+A1.Vp.^2).^0.5,3); +velmagmean=velmagmean(end:-1:1,:); +figure;imagesc(Xp(:),Yp(:),velmagmean); + +% error_prana=[A1.err_up;A1.err_vp]; +% error_davis=[A1.err_ud;A1.err_vd]; +% UMC=[A1.UMCx;A1.UMCy]; +% UIM=[A1.UIMx;A1.UIMy]; +% UCS=[A1.UCSx;A1.UCSy]; +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +[yig,xig]=meshgrid(linspace(12,500,62),linspace(12,1524,190)); + +locx=192/2;% 768 +locx1=locx-1;%764 + +EP=squeeze(error_prana(end:-1:1,locx,:)); +ED=squeeze(error_davis(end:-1:1,locx1,:)); +MC=squeeze(UMC(end:-1:1,locx,:)); +IM=squeeze(UIM(end:-1:1,locx,:)); +CS=squeeze(UCS(end:-1:1,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:62 + k=1; + l=1; + for j=1:100 + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(8:8:503)'; +yax1=(12:8:500)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); + +xax=A1.Xp(1,:)/max(A1.Xp(1,:)); +xpoint=xax(locx)*ones(length(yax),1); +hold on;plot(xpoint,yax,'k--');hold off; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +% plot(RMSerrorprana,yax,'k-'); +% plot(RMSerrordavis,yax1,'k--'); +% plot(RMSUMC,yax,'r-'); +% plot(RMSUIM,yax,'c-'); +% plot(RMSUCS,yax1,'m-'); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +axis([0 1 0 0.2]); + +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'2003B_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'2003B_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS','xpoint','yax','velmagmean','Xp','Yp'); +end + +%% 2005B +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load([basedir,casename,'\matfiles\PivChal05B.mat']); +% Upmean=mean(A1.Up,3); +% figure;imagesc(A1.Xp(:),A1.Yp(:),Upmean); +Xp=A1.Xp/max(A1.Xp(:)); +Yp=A1.Yp/max(A1.Yp(:)); +velmagmean=mean((A1.Up(:,:,1:5).^2+A1.Vp(:,:,1:5).^2).^0.5,3); +% velmagmean=velmagmean(end:-1:1,:); +figure;imagesc(Xp(:),Yp(:),velmagmean(:,:)); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=44;% 712 +locx1=44;%712 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:41 + k=1; + l=1; + for j=1:5 + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(24:16:664)'; +yax1=(24:16:664)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); +xax=A1.Xp(1,:)/max(A1.Xp(1,:)); +xpoint=xax(locx)*ones(length(yax),1); +hold on;plot(xpoint,yax,'k--');hold off; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'2005B_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'2005B_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS','xpoint','yax','velmagmean','Xp','Yp'); +end + +%% stagnation flow +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load([basedir,casename,'\matfiles\stagnation_flow.mat']); +% Upmean=mean(A1.Up,3); +% Vpmean=mean(A1.Vp,3); +% figure;imagesc(A1.Xp(:),A1.Yp(:),Upmean); +% figure;quiver(A1.Xp,A1.Yp,Upmean,Vpmean,5); +% figure;quiver(A1.Xp,A1.Yp,Upmean,Vpmean,5); +Xp=A1.Xp/max(A1.Xp(:)); +Yp=A1.Yp/max(A1.Yp(:)); +% Yp=Yp(end:-1:1,:); +velmagmean=mean((A1.Up.^2+A1.Vp.^2).^0.5,3); +% velmagmean=velmagmean(end:-1:1,:); +figure;imagesc(Xp(:),Yp(:),velmagmean); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=40;% 632 +locx1=40;%632 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:64 + k=1; + l=1; + for j=1:99 + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(8:16:1016)'; +yax1=(8:16:1016)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); +xax=A1.Xp(1,:)/max(A1.Xp(1,:)); +xpoint=xax(locx)*ones(length(yax),1); +hold on;plot(xpoint,yax,'k--');hold off; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +axis([0 1 0 0.2]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'stagnation_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'stagnation_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS','xpoint','yax','velmagmean','Xp','Yp'); +end + +%% Vortex Ring +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load([basedir,casename,'\matfiles\Vortex_Ring.mat']); +% Upmean=mean(A1.Up,3); +% figure;imagesc(A1.Xp(:),A1.Yp(:),Upmean); +Xp=(A1.Xp-min(A1.Xp(:)))/range(A1.Xp(:)); +Yp=(A1.Yp-min(A1.Yp(:)))/range(A1.Yp(:)); +velmagmean=mean((A1.Up.^2+A1.Vp.^2).^0.5,3); +velmagmean=velmagmean(end:-1:1,:); +figure;imagesc(Xp(:),Yp(:),velmagmean); +% figure;imagesc(Xp(:),Yp(:),Upmean); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locy=17;% 632 +locy1=17;%632 + +EP=squeeze(error_prana(locy,:,:)); +ED=squeeze(error_davis(locy,:,:)); +MC=squeeze(UMC(locy,:,:)); +IM=squeeze(UIM(locy,:,:)); +CS=squeeze(UCS(locy,:,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:80 + k=1; + l=1; + for j=1:50 + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(1:1:80)'; +yax1=(1:1:80)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); + +xax=Yp(:,1); +xpoint=xax(locy)*ones(length(yax),1); +hold on;plot(yax,xpoint,'k--');hold off; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +%axis([0 0.2 0 1016]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'vortexring_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'vortexring_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS','xpoint','yax','velmagmean','Xp','Yp'); +end + +%% Jetdata +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load([basedir,casename,'\matfiles\Jetdata.mat']); +% Upmean=mean(A1.Up,3); +% figure;imagesc(A1.Xp(:),A1.Yp(:),Upmean); +Xp=(A1.Xp-min(A1.Xp(:)))/range(A1.Xp(:)); +Yp=-(A1.Yp-254)/(10.2*7.2);%max(A1.Yp(:)); +% Yp=max(A1.Yp(:)); +velmagmean=mean((A1.Up.^2+A1.Vp.^2).^0.5,3); +% velmagmean=velmagmean(end:-1:1,:); +figure;imagesc(Xp(:),Yp(:),velmagmean); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=16;% 632 +locx1=locx;%632 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +% RMSerrorprana=rms(EP,2); +% RMSerrordavis=rms(ED,2); +% RMSUMC=rms(MC,2); +% RMSUIM=rms(IM,2); +% RMSUCS=rms(CS,2); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:31 + k=1; + l=1; + for j=1:495 + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end + +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +% yax=(1:1:31)'; +% yax1=(1:1:31)'; +Yd=A1.Yd(:,1); +yax=-((Yd-254)./(10.2*7.2)); +yax1=yax; +xax=Xp;%A1.Xp(1,:)/max(A1.Xp(1,:)); +xpoint=xax(1,locx)*ones(length(yax),1); +hold on;plot(xpoint,yax,'k--');hold off; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; + axis([-1 1 0 0.2]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'jet_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'jet_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS','xpoint','yax','velmagmean','Xp','Yp'); +end diff --git a/postprocess_codes/RMS_spatial_profiles_2Dmap.m b/postprocess_codes/RMS_spatial_profiles_2Dmap.m new file mode 100644 index 0000000..fd5797e --- /dev/null +++ b/postprocess_codes/RMS_spatial_profiles_2Dmap.m @@ -0,0 +1,420 @@ +%% Plot RMS spatial profiles for error and uncertainty +clear; +casename='WS32'; +outdir=['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\spatial_RMS_profiles_2D\',casename,filesep]; +if ~exist(outdir,'dir'); + mkdir(outdir); +end +savemat=0; +saveplot=0; +% keyboard; +%% 2003B +A1=load(['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\',casename,'\matfiles\PivChal03B.mat']); + +% error_prana=cat(3,A1.err_up,A1.err_vp); +% error_davis=cat(3,A1.err_ud,A1.err_vd); +% UMC=cat(3,A1.UMCx,A1.UMCy); +% UIM=cat(3,A1.UIMx,A1.UIMy); +% UCS=cat(3,A1.UCSx,A1.UCSy); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +[yig,xig]=meshgrid(linspace(12,500,62),linspace(12,1524,190)); + +locx=192/2;% 768 +locx1=locx-1;%764 + +EP=squeeze(error_prana(end:-1:1,locx,:)); +ED=squeeze(error_davis(end:-1:1,locx1,:)); +MC=squeeze(UMC(end:-1:1,locx,:)); +IM=squeeze(UIM(end:-1:1,locx,:)); +CS=squeeze(UCS(end:-1:1,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:62 + k=1; + l=1; + for j=1:size(error_prana,3) + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(8:8:503)'; +yax1=(12:8:500)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +% plot(RMSerrorprana,yax,'k-'); +% plot(RMSerrordavis,yax1,'k--'); +% plot(RMSUMC,yax,'r-'); +% plot(RMSUIM,yax,'c-'); +% plot(RMSUCS,yax1,'m-'); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +axis([0 1 0 0.2]); + +error_prana(abs(error_prana)>thresh)=nan; +rmserrprana=nanmean(error_prana,3); +figure;imagesc(rmserrprana); + +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'2003B_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'2003B_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS'); +end + +%% 2005B +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load(['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\',casename,'\matfiles\PivChal05B.mat']); + +% error_prana=cat(3,A1.err_up,A1.err_vp); +% error_davis=cat(3,A1.err_ud,A1.err_vd); +% UMC=cat(3,A1.UMCx,A1.UMCy); +% UIM=cat(3,A1.UIMx,A1.UIMy); +% UCS=cat(3,A1.UCSx,A1.UCSy); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=44;% 712 +locx1=44;%712 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:41 + k=1; + l=1; + for j=1:size(error_prana,3) + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(24:16:664)'; +yax1=(24:16:664)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'2005B_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'2005B_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS'); +end + +%% stagnation flow +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load(['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\',casename,'\matfiles\stagnation_flow.mat']); + +% error_prana=cat(3,A1.err_up,A1.err_vp); +% error_davis=cat(3,A1.err_ud,A1.err_vd); +% UMC=cat(3,A1.UMCx,A1.UMCy); +% UIM=cat(3,A1.UIMx,A1.UIMy); +% UCS=cat(3,A1.UCSx,A1.UCSy); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=40;% 632 +locx1=40;%632 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:64 + k=1; + l=1; + for j=1:size(error_prana,3) + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(8:16:1016)'; +yax1=(8:16:1016)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +axis([0 1 0 0.2]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'stagnation_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'stagnation_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS'); +end + +%% Vortex Ring +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load(['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\',casename,'\matfiles\Vortex_Ring.mat']); + +% error_prana=cat(3,A1.err_up,A1.err_vp); +% error_davis=cat(3,A1.err_ud,A1.err_vd); +% UMC=cat(3,A1.UMCx,A1.UMCy); +% UIM=cat(3,A1.UIMx,A1.UIMy); +% UCS=cat(3,A1.UCSx,A1.UCSy); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locy=17;% 632 +locy1=17;%632 + +EP=squeeze(error_prana(locy,:,:)); +ED=squeeze(error_davis(locy,:,:)); +MC=squeeze(UMC(locy,:,:)); +IM=squeeze(UIM(locy,:,:)); +CS=squeeze(UCS(locy,:,:)); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:80 + k=1; + l=1; + for j=1:size(error_prana,3) + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +yax=(1:1:80)'; +yax1=(1:1:80)'; +yax=yax/max(yax(:)); +yax1=yax1/max(yax1(:)); + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; +%axis([0 0.2 0 1016]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'vortexring_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'vortexring_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS'); +end + +%% Jetdata +clear RMSerrorprana RMSerrordavis RMSUMC RMSUIM RMSUCS; +clear tempstoreEP tempstoreMC tempstoreIM tempstoreED tempstoreCS; +A1=load(['Z:\Planar_Uncertainty_work\Results\final_plots\03_21_2017\',casename,'\matfiles\Jetdata.mat']); + +% error_prana=cat(3,A1.err_up,A1.err_vp); +% error_davis=cat(3,A1.err_ud,A1.err_vd); +% UMC=cat(3,A1.UMCx,A1.UMCy); +% UIM=cat(3,A1.UIMx,A1.UIMy); +% UCS=cat(3,A1.UCSx,A1.UCSy); + +error_prana=A1.err_up; +error_davis=A1.err_ud; +UMC=A1.UMCx; +UIM=A1.UIMx; +UCS=A1.UCSx; + +locx=16;% 632 +locx1=locx;%632 + +EP=squeeze(error_prana(:,locx,:)); +ED=squeeze(error_davis(:,locx1,:)); +MC=squeeze(UMC(:,locx,:)); +IM=squeeze(UIM(:,locx,:)); +CS=squeeze(UCS(:,locx1,:)); + +% RMSerrorprana=rms(EP,2); +% RMSerrordavis=rms(ED,2); +% RMSUMC=rms(MC,2); +% RMSUIM=rms(IM,2); +% RMSUCS=rms(CS,2); + +thresh=1; +[indx1,indy1]=find(abs(EP)>thresh); +[indx2,indy2]=find(abs(ED)>thresh); + +for i=1:31 + k=1; + l=1; + for j=1:size(error_prana,3) + if ismember(i,indx1) && ismember(j,indy1) + continue; + else + tempstoreEP(k)=EP(i,j); + tempstoreMC(k)=MC(i,j); + tempstoreIM(k)=IM(i,j); + k=k+1; + end + if ismember(i,indx2) && ismember(j,indy2) + continue; + else + tempstoreED(l)=ED(i,j); + tempstoreCS(l)=CS(i,j); + l=l+1; + end + end + RMSerrorprana(i)=rms(tempstoreEP); + RMSerrordavis(i)=rms(tempstoreED); + RMSUMC(i)=rms(tempstoreMC); + RMSUIM(i)=rms(tempstoreIM); + RMSUCS(i)=rms(tempstoreCS); +end + +% yax=[8:8:504 8:8:504]'; +% yax1=[8:8:502 8:8:502]'; +% yax=(1:1:31)'; +% yax1=(1:1:31)'; +Yd=A1.Yd(:,1); +yax=-((Yd-254)./(10.2*7.2)); +yax1=yax; + +lw=2; +figure;hold on;set(gcf,'DefaultLineLineWidth',lw); +plot(yax,RMSerrorprana,'k-'); +plot(yax1,RMSerrordavis,'k--'); +plot(yax,RMSUMC,'r-'); +plot(yax,RMSUIM,'c-'); +plot(yax1,RMSUCS,'m-'); +hold off; + axis([-1 1 0 0.2]) +if saveplot==1 + print(gcf,'-dpng',fullfile(outdir,'jet_spatial_profile.png'),'-r360'); +end +if savemat==1 + save(fullfile(outdir,'jet_spatial_profilemat.mat'),'yax','yax1','RMSerrorprana','RMSerrordavis','RMSUMC','RMSUIM','RMSUCS'); +end diff --git a/postprocess_codes/compute_error_uncertainty.m b/postprocess_codes/compute_error_uncertainty.m new file mode 100644 index 0000000..94c100e --- /dev/null +++ b/postprocess_codes/compute_error_uncertainty.m @@ -0,0 +1,783 @@ +function compute_error_uncertainty(casename,inputdir,outputdir,savemat) +%diameter index +% dn=6;%autod +dn=26;%crossd +% addpath Y:\Projects\TurbulentFlame\analysis\src\readimx_v2.0_win64\readimx; +% addpath Y:\Projects\turbulent_flame\analysis\src\readimx_v2.0_win64\readimx +addpath Z:\Planar_Uncertainty_work\codes\readimx_v2.0_win64\readimx; +if strcmp(casename{1},'PivChal03B') + fprintf('\n PIV Challenge 2003B case uncertainty calculation \n'); + %% Load True Solution + truesol=load(fullfile(inputdir.Pivchal03B.truesoldir,[inputdir.Pivchal03B.truesolbase,'.mat'])); + + %% Set Parameters + N=100; + fstart=1; +% foffset=0; + sxp=63;syp=191; + sxd=62;syd=190; + Up=zeros(sxp,syp,N); + Vp=zeros(sxp,syp,N); + Ud=zeros(sxd,syd,N); + Vd=zeros(sxd,syd,N); + + Utrue=zeros(sxp,syp,N); + Vtrue=zeros(sxp,syp,N); + Utrued=zeros(sxd,syd,N); + Vtrued=zeros(sxd,syd,N); + + err_up=zeros(sxp,syp,N); + err_vp=zeros(sxp,syp,N); + err_ud=zeros(sxd,syd,N); + err_vd=zeros(sxd,syd,N); + + UMCx=zeros(sxp,syp,N); + UMCy=zeros(sxp,syp,N); + UIMx=zeros(sxp,syp,N); + UIMy=zeros(sxp,syp,N); + UCSx=zeros(sxd,syd,N); + UCSy=zeros(sxd,syd,N); + Udiffstore=zeros(sxp,syp,N); + Vdiffstore=zeros(sxp,syp,N); + %% Loop through N fields + for i=fstart:N + + %% Load prana solution + pranasol=load(fullfile(inputdir.Pivchal03B.pranasol,[inputdir.Pivchal03B.pranabaseMC, num2str(2*i-1,'%03.0f'),'.mat'])); + imsol=load(fullfile(inputdir.Pivchal03B.IMdir,[inputdir.Pivchal03B.pranabaseIM, num2str(2*i-1,'%03.0f'),'.mat'])); +% keyboard; + Up(:,:,i)=pranasol.U(end:-1:1,:,1); + Vp(:,:,i)=pranasol.V(end:-1:1,:,1); + if i==fstart + Xp=pranasol.X; + Yp=pranasol.Y; + end + %% Load Davis Solution + davissol=readimx(fullfile(inputdir.Pivchal03B.davisol,[inputdir.Pivchal03B.davisbase,num2str(i,'%05.0f'),'.vc7'])); + Utemp=cell2mat(davissol.Frames{1,1}.Components{1,1}.Planes)'; + Vtemp=cell2mat(davissol.Frames{1,1}.Components{2,1}.Planes)'; + Ud(:,:,i)=Utemp(2:end-1,2:end-1); + Vd(:,:,i)=-Vtemp(2:end-1,2:end-1); + %Davis X and Y grid + if i==fstart + [D] = create2DVec(davissol.Frames{1,1}); + Xd=D.X(2:end-1,2:end-1); + Yd=D.Y(2:end-1,2:end-1); + clear D; + Xt=truesol.x; + Yt=truesol.y(:,end:-1:1); + + [yig,xig]=meshgrid(linspace(12,500,62),linspace(12,1524,190)); + + end + + %% Get True Solution + Utrue(:,:,i)=truesol.u(:,:,i)'; + Vtrue(:,:,i)=truesol.v(:,:,i)'; + %Interpolated true solution onto Davis Grid + Utrued(:,:,i)=interp2(Xt',Yt',Utrue(:,:,i),xig',yig','linear',0); + Vtrued(:,:,i)=interp2(Xt',Yt',Vtrue(:,:,i),xig',yig','linear',0); + %% Calculate Error + err_up(:,:,i)=(Up(:,:,i)-Utrue(:,:,i)); + err_vp(:,:,i)=(Vp(:,:,i)-Vtrue(:,:,i)); + err_ud(:,:,i)=(Ud(:,:,i)-Utrued(:,:,i)); + err_vd(:,:,i)=(Vd(:,:,i)-Vtrued(:,:,i)); + + %% Get MC uncertainty estimates + Ixx=pranasol.ixx(end:-1:1,:); + Iyy=pranasol.iyy(end:-1:1,:); + scaling=pranasol.mi(end:-1:1,:); + biasx=reshape(pranasol.uncertainty(:,15),sxp,syp); + biasy=reshape(pranasol.uncertainty(:,16),sxp,syp); + Autod=reshape(pranasol.uncertainty(:,dn),sxp,syp); + biasx=biasx(end:-1:1,:); + biasy=biasy(end:-1:1,:); + Autod=Autod(end:-1:1,:); + + %Gradient correction + Udiff=socdiff(Up(:,:,i),8,1); + Vdiff=socdiff(Vp(:,:,i),8,2); +% [dudx ]=gradient_compact_rich(N,dx,dir) +% Udiff=gradient_compact_rich(Up(:,:,i),8,1); +% Vdiff=gradient_compact_rich(Vp(:,:,i),8,2); + Udiffstore(:,:,i)=Udiff; + Vdiffstore(:,:,i)=Vdiff; + + Ixxt= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 )); + Iyyt= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 )); + %scaling and bias correction + UMCx(:,:,i)=sqrt(biasx.^2+(Ixxt.^2)./scaling); + UMCy(:,:,i)=sqrt(biasy.^2+(Iyyt.^2)./scaling); + + + %% Get IM uncertainty estimates + UIMx(:,:,i)=imsol.imx(end:-1:1,:); + UIMy(:,:,i)=imsol.imy(end:-1:1,:); + + %% Get CS uncertainty estimates + Unrxtemp=cell2mat(davissol.Frames{1,1}.Components{22,1}.Planes)'; + Unrytemp=cell2mat(davissol.Frames{1,1}.Components{23,1}.Planes)'; + Unbxtemp=cell2mat(davissol.Frames{1,1}.Components{25,1}.Planes)'; + Unbytemp=cell2mat(davissol.Frames{1,1}.Components{26,1}.Planes)'; + + Unrx=Unrxtemp(2:end-1,2:end-1); + Unry=Unrytemp(2:end-1,2:end-1); + Unbx=Unbxtemp(2:end-1,2:end-1); + Unby=Unbytemp(2:end-1,2:end-1); + + UCSx(:,:,i)=sqrt(Unrx.^2+Unbx.^2); + UCSy(:,:,i)=sqrt(Unry.^2+Unby.^2); + + + end + Xd=xig'; + Yd=yig'; + %% Save error and uncertainty values + if savemat==1 + save(fullfile(outputdir{1},'PivChal03B.mat'),'err_up','err_vp','err_ud','err_vd','UMCx','UMCy','UIMx','UIMy','UCSx','UCSy','Up','Vp','Ud','Vd','Xp','Yp','Xt','Yt','Xd','Yd','Udiff','Vdiff'); + end + clearvars -except casename inputdir outputdir savemat dn; + fprintf('\nDone...\n'); + +end +%% +if strcmp(casename{2},'PivChal05B') + fprintf('\n PIV Challenge 2005B case uncertainty calculation \n'); + %% Load True Solution + truesol=load(fullfile(inputdir.Pivchal05B.truesoldir,[inputdir.Pivchal05B.truesolbase,'.mat'])); + + %% Set Parameters + N=6; + fstart=1; +% foffset=0; + f=[9 29 49 69 89 109]; + + sxp=41;syp=88; + sxd=41;syd=88; + Up=zeros(sxp,syp,N); + Vp=zeros(sxp,syp,N); + Ud=zeros(sxd,syd,N); + Vd=zeros(sxd,syd,N); + + Utrue=zeros(sxp,syp,N); + Vtrue=zeros(sxp,syp,N); +% Utrued=zeros(sxd,syd,N); +% Vtrued=zeros(sxd,syd,N); + + err_up=zeros(sxp,syp,N); + err_vp=zeros(sxp,syp,N); + err_ud=zeros(sxd,syd,N); + err_vd=zeros(sxd,syd,N); + + UMCx=zeros(sxp,syp,N); + UMCy=zeros(sxp,syp,N); + UIMx=zeros(sxp,syp,N); + UIMy=zeros(sxp,syp,N); + UCSx=zeros(sxd,syd,N); + UCSy=zeros(sxd,syd,N); + Udiffstore=zeros(sxp,syp,N); + Vdiffstore=zeros(sxp,syp,N); + + %% Loop through N fields + for i=fstart:N + + %% Load prana solution + pranasol=load(fullfile(inputdir.Pivchal05B.pranasol,[inputdir.Pivchal05B.pranabaseMC, num2str(f(i),'%03.0f'),'.mat'])); + imsol=load(fullfile(inputdir.Pivchal05B.IMdir,[inputdir.Pivchal05B.pranabaseIM, num2str(f(i),'%03.0f'),'.mat'])); +% keyboard; + Up(:,:,i)=pranasol.U./2; + Vp(:,:,i)=pranasol.V./2; + if i==fstart + Xp=pranasol.X; + Yp=pranasol.Y; + end + %% Load Davis Solution + davissol=readimx(fullfile(inputdir.Pivchal05B.davisol,[inputdir.Pivchal05B.davisbase,num2str(i,'%05.0f'),'.vc7'])); +% keyboard; + Utemp=cell2mat(davissol.Frames{1,1}.Components{1,1}.Planes)'; + Vtemp=cell2mat(davissol.Frames{1,1}.Components{2,1}.Planes)'; + Ud(:,:,i)=Utemp(end-1:-1:2,2:end-1)./2; + Vd(:,:,i)=-Vtemp(end-1:-1:2,2:end-1)./2; + %Davis X and Y grid + if i==fstart + [D] = create2DVec(davissol.Frames{1,1}); + Xd=D.X(end-1:-1:2,2:end-1); + Yd=D.Y(end-1:-1:2,2:end-1); + clear D; + Xt=truesol.Xt; + Yt=truesol.Yt; + end + + %% Get True Solution + Utrue(:,:,i)=truesol.Ut(:,:,i); + Vtrue(:,:,i)=truesol.Vt(:,:,i); + + %% Calculate Error + err_up(:,:,i)=(Up(:,:,i)-Utrue(:,:,i)); + err_vp(:,:,i)=(Vp(:,:,i)-Vtrue(:,:,i)); + err_ud(:,:,i)=(Ud(:,:,i)-Utrue(:,:,i)); + err_vd(:,:,i)=(Vd(:,:,i)-Vtrue(:,:,i)); + + %% Get MC uncertainty estimates + Ixx=pranasol.ixx./2; + Iyy=pranasol.iyy./2; + scaling=pranasol.mi; + biasx=reshape(pranasol.uncertainty(:,15),sxp,syp)./2; + biasy=reshape(pranasol.uncertainty(:,16),sxp,syp)./2; + Autod=reshape(pranasol.uncertainty(:,dn),sxp,syp); + + %Gradient correction + Udiff=socdiff(Up(:,:,i),8,1); + Vdiff=socdiff(Vp(:,:,i),8,2); +% Udiff=gradient_compact_rich(Up(:,:,i),8,1); +% Vdiff=gradient_compact_rich(Vp(:,:,i),8,2); + Udiffstore(:,:,i)=Udiff; + Vdiffstore(:,:,i)=Vdiff; + + Ixxt= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 )); + Iyyt= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 )); + %scaling and bias correction + UMCx(:,:,i)=sqrt(biasx.^2+(Ixxt.^2)./scaling); + UMCy(:,:,i)=sqrt(biasy.^2+(Iyyt.^2)./scaling); + + + %% Get IM uncertainty estimates + UIMx(:,:,i)=imsol.imx./2; + UIMy(:,:,i)=imsol.imy./2; + + %% Get CS uncertainty estimates + Unrxtemp=cell2mat(davissol.Frames{1,1}.Components{22,1}.Planes)'; + Unrytemp=cell2mat(davissol.Frames{1,1}.Components{23,1}.Planes)'; + Unbxtemp=cell2mat(davissol.Frames{1,1}.Components{25,1}.Planes)'; + Unbytemp=cell2mat(davissol.Frames{1,1}.Components{26,1}.Planes)'; + + Unrx=Unrxtemp(end-1:-1:2,2:end-1)./2; + Unry=Unrytemp(end-1:-1:2,2:end-1)./2; + Unbx=Unbxtemp(end-1:-1:2,2:end-1)./2; + Unby=Unbytemp(end-1:-1:2,2:end-1)./2; + + UCSx(:,:,i)=sqrt(Unrx.^2+Unbx.^2); + UCSy(:,:,i)=sqrt(Unry.^2+Unby.^2); + + + end + %% Save error and uncertainty values + if savemat==1 + save(fullfile(outputdir{1},'PivChal05B.mat'),'err_up','err_vp','err_ud','err_vd','UMCx','UMCy','UIMx','UIMy','UCSx','UCSy','Up','Vp','Ud','Vd','Xp','Yp','Xt','Yt','Xd','Yd','Udiff','Vdiff'); + end + clearvars -except casename inputdir outputdir savemat dn; + fprintf('\nDone...\n'); + +end + +%% +if strcmp(casename{3},'stagnation_flow') + fprintf('\n stagnation_flow case uncertainty calculation \n'); + %% Load True Solution + truesol=load(fullfile(inputdir.stagflow.truesoldir,[inputdir.stagflow.truesolbase,'.mat'])); + + %% Set Parameters + N=99; + fstart=1; + + sxp=64;syp=80; + sxd=64;syd=80; + Up=zeros(sxp,syp,N); + Vp=zeros(sxp,syp,N); + Ud=zeros(sxd,syd,N); + Vd=zeros(sxd,syd,N); + + Utrue=zeros(sxp,syp,N); + Vtrue=zeros(sxp,syp,N); +% Utrued=zeros(sxd,syd,N); +% Vtrued=zeros(sxd,syd,N); + + err_up=zeros(sxp,syp,N); + err_vp=zeros(sxp,syp,N); + err_ud=zeros(sxd,syd,N); + err_vd=zeros(sxd,syd,N); + + UMCx=zeros(sxp,syp,N); + UMCy=zeros(sxp,syp,N); + UIMx=zeros(sxp,syp,N); + UIMy=zeros(sxp,syp,N); + UCSx=zeros(sxd,syd,N); + UCSy=zeros(sxd,syd,N); + Udiffstore=zeros(sxp,syp,N); + Vdiffstore=zeros(sxp,syp,N); + + %% Loop through N fields + for i=fstart:N + + %% Load prana solution + pranasol=load(fullfile(inputdir.stagflow.pranasol,[inputdir.stagflow.pranabaseMC, num2str(i,'%06.0f'),'.mat'])); + imsol=load(fullfile(inputdir.stagflow.IMdir,[inputdir.stagflow.pranabaseIM, num2str(i,'%06.0f'),'.mat'])); +% keyboard; + Up(:,:,i)=pranasol.U; + Vp(:,:,i)=pranasol.V; + if i==fstart + Xp=pranasol.X; + Yp=pranasol.Y; + end + %% Load Davis Solution + davissol=readimx(fullfile(inputdir.stagflow.davisol,[inputdir.stagflow.davisbase,num2str(i,'%05.0f'),'.vc7'])); +% keyboard; + Utemp=cell2mat(davissol.Frames{1,1}.Components{1,1}.Planes)'; + Vtemp=cell2mat(davissol.Frames{1,1}.Components{2,1}.Planes)'; + + Ud(:,:,i)=Utemp(end:-1:1,:); + Vd(:,:,i)=-Vtemp(end:-1:1,:); + %Davis X and Y grid + if i==fstart + [D] = create2DVec(davissol.Frames{1,1}); + Xd=D.X(:,:)'; +% Xd=Xd(:,end:-1:1); + Yd=D.Y(end:-1:1,:)'; + clear D; +% Xt=truesol.Xt; +% Yt=truesol.Yt; + end + + %% Get True Solution + Utrue(:,:,i)=truesol.Ut; + Vtrue(:,:,i)=truesol.Vt; + if i==fstart + unfitx=truesol.unfitx; + unfity=truesol.unfity; + end + %% Calculate Error + err_up(:,:,i)=(Up(:,:,i)-Utrue(:,:,i)); + err_vp(:,:,i)=(Vp(:,:,i)-Vtrue(:,:,i)); + err_ud(:,:,i)=(Ud(:,:,i)-Utrue(:,:,i)); + err_vd(:,:,i)=(Vd(:,:,i)-Vtrue(:,:,i)); + + %% Get MC uncertainty estimates + Ixx=pranasol.ixx; + Iyy=pranasol.iyy; + scaling=pranasol.mi; + biasx=reshape(pranasol.uncertainty(:,15),sxp,syp); + biasy=reshape(pranasol.uncertainty(:,16),sxp,syp); + Autod=reshape(pranasol.uncertainty(:,dn),sxp,syp); + + %Gradient correction + Udiff=socdiff(Up(:,:,i),16,1); + Vdiff=socdiff(Vp(:,:,i),16,2); +% Udiff=gradient_compact_rich(Up(:,:,i),16,1); +% Vdiff=gradient_compact_rich(Vp(:,:,i),16,2); + Udiffstore(:,:,i)=Udiff; + Vdiffstore(:,:,i)=Vdiff; + + Ixxt= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 )); + Iyyt= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 )); + %scaling and bias correction + UMCx(:,:,i)=sqrt(biasx.^2+(Ixxt.^2)./scaling +unfitx.^2); + UMCy(:,:,i)=sqrt(biasy.^2+(Iyyt.^2)./scaling +unfity.^2); + + + %% Get IM uncertainty estimates + UIMx(:,:,i)=sqrt((imsol.imx).^2 +unfitx.^2); + UIMy(:,:,i)=sqrt((imsol.imy).^2 +unfity.^2); + + %% Get CS uncertainty estimates + Unrxtemp=cell2mat(davissol.Frames{1,1}.Components{22,1}.Planes)'; + Unrytemp=cell2mat(davissol.Frames{1,1}.Components{23,1}.Planes)'; + Unbxtemp=cell2mat(davissol.Frames{1,1}.Components{25,1}.Planes)'; + Unbytemp=cell2mat(davissol.Frames{1,1}.Components{26,1}.Planes)'; + + Unrx=Unrxtemp(end:-1:1,:); + Unry=Unrytemp(end:-1:1,:); + Unbx=Unbxtemp(end:-1:1,:); + Unby=Unbytemp(end:-1:1,:); + + UCSx(:,:,i)=sqrt(Unrx.^2+Unbx.^2 +unfitx.^2); + UCSy(:,:,i)=sqrt(Unry.^2+Unby.^2 +unfity.^2); +% keyboard; + + end + Xt=Xp; + Yt=Yp; + %% Save error and uncertainty values + if savemat==1 + save(fullfile(outputdir{1},'stagnation_flow.mat'),'err_up','err_vp','err_ud','err_vd','UMCx','UMCy','UIMx','UIMy','UCSx','UCSy','Up','Vp','Ud','Vd','Xp','Yp','Xt','Yt','Xd','Yd','Udiff','Vdiff'); + end + clearvars -except casename inputdir outputdir savemat dn; + fprintf('\nDone...\n'); + +end + +%% +if strcmp(casename{4},'Vortex_Ring') + fprintf('\n Vortex_Ring case uncertainty calculation \n'); + %% Load True Solution + truesol=load(fullfile(inputdir.vortex.truesoldir,[inputdir.vortex.truesolbase,'.mat'])); + + %% Set Parameters + N=50; + fstart=1; + foffsetprana=24; + foffsetdavis=1; + + spf=1; + xscale=1/11.9690;%stereojob.datasave.scaling.wil.xscale; + yscale=1/11.9653;%stereojob.datasave.scaling.wil.yscale; + xorigin=576.8623; + yorigin=395.8564; + + %true solution grid + [Xt,Yt]=meshgrid(-34:0.45:20,-27:0.45:27); + %mapping true solution grid to pixel grid using magnification and + %origin shift + Xt1=Xt./xscale+xorigin; + Yt1=Yt./yscale+yorigin; + % get limits of true solution grid in pixel coordinates + ax=min(Xt1(:)); + bx=max(Xt1(:)); + ay=min(Yt1(:)); + by=max(Yt1(:)); + + %Load 1 prana sol + soltemp=load(fullfile(inputdir.vortex.pranasol,[inputdir.vortex.pranabaseMC, num2str(1+foffsetprana,'%05.0f'),'.mat'])); + %prana solution grid + Xsol1=soltemp.X; + Ysol1=soltemp.Y; + %Load 1 Davis sol + davissoltemp=readimx(fullfile(inputdir.vortex.davisol,[inputdir.vortex.davisbase,num2str(1+foffsetdavis,'%05.0f'),'.vc7'])); + %davis solution grid + D=create2DVec(davissoltemp.Frames{1,1}); + Xd1=D.X(:,:)'; + Yd1=D.Y(end:-1:1,:)'; + clear D soltemp davissoltemp; + % indidces corresponding to region of interest + + xvec=Xsol1(1,:); + yvec=Ysol1(:,1); + xvec1=Xd1(1,:); + yvec1=Yd1(:,1); + + %Get indices which will be considered from the prana solution + xind=find(xvec>ax & xvecay & yvecax & xvec1ay & yvec1 3 + error('function only defined up to 3D matrices') +end + +if isempty(dx) + dx = 1; +end + + +N = permute(N, [dir, 1:dir-1, dir+1:ndim_N]); +size_N = size(N); %find again for permuted matrix +ndim_N = ndims(N); + + +NI = size_N(1); +NJ = size_N(2); +if ndim_N == 3 + NK = size_N(3); +else + NK = 1; +end +% DN = zeros(NI,NJ,NK); +ADN = zeros(NI,NJ,NK); + + +A=[1239 272 1036 -69 0]; +k=[1 2 4 8]; + +%% Parameters for Boundary formulation for the First +%% Derivative (Lele et al 1992) +alpha=1; +d1=0; +a1=-(3+alpha+2*d1)/2; +b1=2+3*d1; +c1=-(1-alpha+6*d1)/2; + + +for j=1:NJ + current_count=1; + for m=1:4; + c1=0; + while c1<=k + N_current(:,1)=N(k(m)+c1:k(m):size(N,1),j); + + a(:,1)=1/4*ones(1,size(N_current,1)); + a(1)=0; + a(size(N_current,1))=0; + b(:,1)=ones(1,size(N_current,1)); + b(1)=1; + b(size(N_current,1))=1; + c(:,1)=1/4*ones(1,size(N_current,1)); + c(1)=0; + c(size(N_current,1))=0; + for l=2:size(N_current,1)-1 + d(l,1)=1.5*(N_current(l+1,1)-N_current(l-1,1))/(2*k(m)*dx); + end + % d(1,1)=1/(k(m)*dx)*(a1*N_current(1,1)+b1*N_current(2,1)+c1*N_current(3,1)+d1*N_current(4,1)); + % d(size(N_current,1),1)=1/(k(m)*dx)*(a1*N_current(size(N_current,1),1)+b1*N_current(size(N_current,1)-1,1)+c1*N_current(size(N_current,1)-2,1)-d1*N_current(size(N_current,1)-3,1)); + % % d(1,1)=0; + % d(size(N_current,1),1)=0; + + d(1,1)=1/(k(m)*dx)*(N_current(2,1)-N_current(1,1)); + d(size(N_current,1),1)=1/(k(m)*dx)*(N_current(size(N_current,1),1)-N_current(size(N_current,1)-1,1)); + % + + [y]=tridiagSolve(a,b,c,d); %% Solve tridiagonal Matrix + ADN(k(m)+c1:k(m):size(N,1),j)=A(m+1)*y'+ADN(k(m)+c1:k(m):size(N,1),j); %% Calculate current sum of (A*U') + % pause + clear a b c d y + + clear N_current + c1=c1+1; + end + end + DN=1/A(1)*ADN; +end +dudx = ipermute(DN, [dir, 1:dir-1, dir+1:ndim_N]); diff --git a/postprocess_codes/gradient_compactrich.m b/postprocess_codes/gradient_compactrich.m new file mode 100644 index 0000000..119ce0f --- /dev/null +++ b/postprocess_codes/gradient_compactrich.m @@ -0,0 +1,133 @@ +function [dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz] = gradient_compactrich(u,v,w,dx,dy,dz) + +for i=1:size(u,3) + dudx(:,:,i)=subfunction_compactrich(u(:,:,i),dx,2); + dudy(:,:,i)=subfunction_compactrich(u(:,:,i),dy,1); + dvdx(:,:,i)=subfunction_compactrich(v(:,:,i),dx,2); + dvdy(:,:,i)=subfunction_compactrich(v(:,:,i),dy,1); + dwdx(:,:,i)=subfunction_compactrich(w(:,:,i),dx,2); + dwdy(:,:,i)=subfunction_compactrich(w(:,:,i),dy,1); +end +for i=1:size(u,2) + dudz(:,i,:)=subfunction_compactrich(u(:,i,:),dz,3); + dvdz(:,i,:)=subfunction_compactrich(v(:,i,:),dz,3); + dwdz(:,i,:)=subfunction_compactrich(w(:,i,:),dz,3); +end + + function [dudx]=subfunction_compactrich(N,dx,dir) + % This function calculates gradients using the 4th order compact-richardson + % scheme introduced by A. Etebari and P. Vlachos in "Improvements on + % the accuracy of derivative estimation from DPIV velocity measurements" + % Experiments in Fluids (2005) + + size_N = size(N); + + if isempty(dir) + if size_N(1)~=1 + dir = 1; + elseif size_N(2)~=1 + dir = 2; + elseif size_N(3)~=1 + dir = 3; + else + dir = 1; + end + end + + ndim_N = ndims(N); + if ndim_N > 3 + error('function only defined up to 3D matrices') + end + + if isempty(dx) + dx = 1; + end + + + N = permute(N, [dir, 1:dir-1, dir+1:ndim_N]); + size_N = size(N); %find again for permuted matrix + ndim_N = ndims(N); + + + NI = size_N(1); + NJ = size_N(2); + if ndim_N == 3 + NK = size_N(3); + else + NK = 1; + end + % DN = zeros(NI,NJ,NK); + ADN = zeros(NI,NJ,NK); + + + A=[1239 272 1036 -69 0]; + k=[1 2 4 8]; + + % Parameters for Boundary formulation for the First + % Derivative (Lele et al 1992) + alpha=1; + d1=0; + a1=-(3+alpha+2*d1)/2; + b1=2+3*d1; + c1=-(1-alpha+6*d1)/2; + + % alpha=2; + % d1=0; + % a1=-(11+2*alpha)/6; + % b1=(6-alpha)/2; + % c1=-(2*alpha-3)/2; + % for i=1:NI; + + for j=1:NJ + % for j=110; + % for j=1; + % current_count=1; + for m=1:4; + + % while current_count<=k; + c1=0; + while c1<=k + N_current(:,1)=N(k(m)+c1:k(m):size(N,1),j); + + a(:,1)=1/4*ones(1,size(N_current,1)); + a(1)=0; + a(size(N_current,1))=0; + b(:,1)=ones(1,size(N_current,1)); + b(1)=1; + b(size(N_current,1))=1; + c(:,1)=1/4*ones(1,size(N_current,1)); + c(1)=0; + c(size(N_current,1))=0; + for l=2:size(N_current,1)-1 + d(l,1)=1.5*(N_current(l+1,1)-N_current(l-1,1))/(2*k(m)*dx); + end + % d(1,1)=1/(k(m)*dx)*(a1*N_current(1,1)+b1*N_current(2,1)+c1*N_current(3,1)+d1*N_current(4,1)); + % d(size(N_current,1),1)=1/(k(m)*dx)*(a1*N_current(size(N_current,1),1)+b1*N_current(size(N_current,1)-1,1)+c1*N_current(size(N_current,1)-2,1)-d1*N_current(size(N_current,1)-3,1)); + % % d(1,1)=0; + % d(size(N_current,1),1)=0; + + d(1,1)=1/(k(m)*dx)*(N_current(2,1)-N_current(1,1)); + d(size(N_current,1),1)=1/(k(m)*dx)*(N_current(size(N_current,1),1)-N_current(size(N_current,1)-1,1)); + % + + [y]=tridiagSolve(a,b,c, d); %% Solve tridiagonal Matrix + ADN(k(m)+c1:k(m):size(N,1),j)=A(m+1)*y'+ADN(k(m)+c1:k(m):size(N,1),j); %% Calculate current sum of (A*U') + % pause + clear a b c d y + + + clear N_current + c1=c1+1; + end + % current_count=current_count+1; + % end + % % + % % + % % + end + DN=1/A(1)*ADN; + end + % end + dudx = ipermute(DN, [dir, 1:dir-1, dir+1:ndim_N]); + end +end \ No newline at end of file diff --git a/postprocess_codes/gradient_dudxdvdy_compact_rich.m b/postprocess_codes/gradient_dudxdvdy_compact_rich.m new file mode 100644 index 0000000..0310f37 --- /dev/null +++ b/postprocess_codes/gradient_dudxdvdy_compact_rich.m @@ -0,0 +1,129 @@ +function [dudx ]=gradient_dudxdvdy_compact_rich(N,dx,dir) + +% This function calculates gradients using the 4th order compact-richardson +% scheme introduced by A. Etebari and P. Vlachos in "Improvements on +% the accuracy of derivative estimation from DPIV velocity measurements" +% Experiments in Fluids (2005) + +% % clc +% clear N_current a b c d A row col DN +% % % +% % % N=V1_recon(:,:,140); +% % % N1=N; +% N=Vs_d; +% dir=2; +% dx=0.000248 +% % % +% % % dir=1; +% % % % dx=0.000248; %m +% % % dx=1; %m + +size_N = size(N); + +if isempty(dir) + if size_N(1)~=1 + dir = 1; + elseif size_N(2)~=1 + dir = 2; + elseif size_N(3)~=1 + dir = 3; + else + dir = 1; + end +end + +ndim_N = ndims(N); +if ndim_N > 3 + error('function only defined up to 3D matrices') +end + +if isempty(dx) + dx = 1; +end + + +N = permute(N, [dir, 1:dir-1, dir+1:ndim_N]); +size_N = size(N); %find again for permuted matrix +ndim_N = ndims(N); + + +NI = size_N(1); +NJ = size_N(2); +if ndim_N == 3 + NK = size_N(3); +else + NK = 1; +end +% DN = zeros(NI,NJ,NK); +ADN = zeros(NI,NJ,NK); + + +A=[1239 272 1036 -69 0]; +k=[1 2 4 8]; + + %% Parameters for Boundary formulation for the First + %% Derivative (Lele et al 1992) + alpha=1; + d1=0; + a1=-(3+alpha+2*d1)/2; + b1=2+3*d1; + c1=-(1-alpha+6*d1)/2; + +% alpha=2; +% d1=0; +% a1=-(11+2*alpha)/6; +% b1=(6-alpha)/2; +% c1=-(2*alpha-3)/2; +% for i=1:NI; + + for j=1:NJ +% for j=110; +% for j=1; + current_count=1; + for m=1:4; + +% while current_count<=k; + c1=0; + while c1<=k + N_current(:,1)=N(k(m)+c1:k(m):size(N,1),j); + + a(:,1)=1/4*ones(1,size(N_current,1)); + a(1)=0; + a(size(N_current,1))=0; + b(:,1)=ones(1,size(N_current,1)); + b(1)=1; + b(size(N_current,1))=1; + c(:,1)=1/4*ones(1,size(N_current,1)); + c(1)=0; + c(size(N_current,1))=0; + for l=2:size(N_current,1)-1 + d(l,1)=1.5*(N_current(l+1,1)-N_current(l-1,1))/(2*k(m)*dx); + end +% d(1,1)=1/(k(m)*dx)*(a1*N_current(1,1)+b1*N_current(2,1)+c1*N_current(3,1)+d1*N_current(4,1)); +% d(size(N_current,1),1)=1/(k(m)*dx)*(a1*N_current(size(N_current,1),1)+b1*N_current(size(N_current,1)-1,1)+c1*N_current(size(N_current,1)-2,1)-d1*N_current(size(N_current,1)-3,1)); +% % d(1,1)=0; +% d(size(N_current,1),1)=0; + + d(1,1)=1/(k(m)*dx)*(N_current(2,1)-N_current(1,1)); + d(size(N_current,1),1)=1/(k(m)*dx)*(N_current(size(N_current,1),1)-N_current(size(N_current,1)-1,1)); +% + + [y]=tridiagSolve(a,b,c, d); %% Solve tridiagonal Matrix + ADN(k(m)+c1:k(m):size(N,1),j)=A(m+1)*y'+ADN(k(m)+c1:k(m):size(N,1),j); %% Calculate current sum of (A*U') +% pause + clear a b c d y + + + clear N_current + c1=c1+1; + end +% current_count=current_count+1; +% end +% % +% % +% % + end + DN=1/A(1)*ADN; + end +% end +dudx = ipermute(DN, [dir, 1:dir-1, dir+1:ndim_N]); diff --git a/postprocess_codes/in_out_directory_set.m b/postprocess_codes/in_out_directory_set.m new file mode 100644 index 0000000..caa7e15 --- /dev/null +++ b/postprocess_codes/in_out_directory_set.m @@ -0,0 +1,79 @@ +function [inputdir,outputdir]=in_out_directory_set(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +% runname='10_14_2015_deform_whittaker_blackman'; + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(4x64x64_87%ov)_01'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname,'2003B','ws32IM'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(4x64x64_75%ov)'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname,'2005B','ws32IM'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(2x32x32_50%ov)'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname,'stagnation_flow','ws32IM'); +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(4x64x64_87%ov)'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname,'Vortex_Ring','ws32IM'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; + +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303_y_188_frame_03_501','TR_PIV_MP(3x16x16_75%ov)_01'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname,'Jetdata','ws32IM'); +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/in_out_directory_set_03_09_2017.m b/postprocess_codes/in_out_directory_set_03_09_2017.m new file mode 100644 index 0000000..15157d5 --- /dev/null +++ b/postprocess_codes/in_out_directory_set_03_09_2017.m @@ -0,0 +1,80 @@ +function [inputdir,outputdir]=in_out_directory_set_03_09_2017(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +runname2='10_22_2015_different_processing_windows';% This for image matching old processing directory + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname,'SOC_gradient'); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure +% runname2='10_22_2015_different_processing_windows'; + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname2,'2003B'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname2,'2005B'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname2,'stagnation_flow'); +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname2,'Vortex_Ring'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; +% +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501','TR_PIV_MP(1x32x32_75%ov)(3x16x16_75%ov)_nopostprocess'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname2,'Jetdata'); +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/in_out_directory_set_03_09_2017_ws64.m b/postprocess_codes/in_out_directory_set_03_09_2017_ws64.m new file mode 100644 index 0000000..2630f88 --- /dev/null +++ b/postprocess_codes/in_out_directory_set_03_09_2017_ws64.m @@ -0,0 +1,79 @@ +function [inputdir,outputdir]=in_out_directory_set_03_09_2017_ws64(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +runname2='10_14_2015_deform_whittaker_blackman';% This for image matching old processing directory + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname,'SOC_gradient'); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(4x64x64_87%ov)_01'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname2,'2003B','ws32IM'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(4x64x64_75%ov)'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname2,'2005B','ws32IM'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(4x64x64_75%ov)'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname,'stagnation_flow');% New IM processing done +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(4x64x64_87%ov)'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname2,'Vortex_Ring','ws32IM'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; + +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501','TR_PIV_MP(4x32x32_87%ov)_nopostprocess'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname,'Jetdata');% New IM processing done +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/in_out_directory_set_07_30_2018_ws1.m b/postprocess_codes/in_out_directory_set_07_30_2018_ws1.m new file mode 100644 index 0000000..6a660ba --- /dev/null +++ b/postprocess_codes/in_out_directory_set_07_30_2018_ws1.m @@ -0,0 +1,79 @@ +function [inputdir,outputdir]=in_out_directory_set_07_30_2018_ws1(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +% runname2='10_14_2015_deform_whittaker_blackman';% This for image matching old processing directory + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname,'SOC_gradient'); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(4x64x64_87%ov)_01'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(4x64x64_75%ov)'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(4x64x64_75%ov)'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname,'stagnation_flow');% New IM processing done +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(4x64x64_87%ov)'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; + +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501','TR_PIV_MP(4x32x32_87%ov)_nopostprocess'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname,'Jetdata');% New IM processing done +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/in_out_directory_set_07_30_2018_ws2.m b/postprocess_codes/in_out_directory_set_07_30_2018_ws2.m new file mode 100644 index 0000000..932bfe9 --- /dev/null +++ b/postprocess_codes/in_out_directory_set_07_30_2018_ws2.m @@ -0,0 +1,80 @@ +function [inputdir,outputdir]=in_out_directory_set_07_30_2018_ws2(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +% runname2='10_22_2015_different_processing_windows';% This for image matching old processing directory + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname,'SOC_gradient'); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure +% runname2='10_22_2015_different_processing_windows'; + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; +% +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501','TR_PIV_MP(1x32x32_75%ov)(3x16x16_75%ov)_nopostprocess'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/in_out_directory_set_10_22_2015.m b/postprocess_codes/in_out_directory_set_10_22_2015.m new file mode 100644 index 0000000..afaa050 --- /dev/null +++ b/postprocess_codes/in_out_directory_set_10_22_2015.m @@ -0,0 +1,79 @@ +function [inputdir,outputdir]=in_out_directory_set_10_22_2015(basedirname,runname) +%Set Input Output Directories + +% The case for which diectories will be created +% runname='10_14_2015_deform_whittaker_blackman'; + +%Creating outputdirectory structure +outputdir_main=fullfile(basedirname,'Results','final_plots',runname); +if ~exist(outputdir_main,'dir') + mkdir(outputdir_main); +end +outputdir{1}=fullfile(outputdir_main,'matfiles'); +outputdir{2}=fullfile(outputdir_main,'figures'); +if ~exist(outputdir{1},'dir') + mkdir(outputdir{1}); +end +if ~exist(outputdir{2},'dir') + mkdir(outputdir{2}); +end + +% Get input directory structure + +%2003B +inputdir.Pivchal03B.pranasol=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2003B','2003bimages','CreateMultiframe','PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.Pivchal03B.MCdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.IMdir=fullfile(basedirname,'Results',runname,'2003B'); +inputdir.Pivchal03B.pranabaseMC='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.pranabaseIM='PIV_scc_finalchal03_pass4_'; +inputdir.Pivchal03B.davisbase='B'; +inputdir.Pivchal03B.truesoldir=fullfile(basedirname,'Results','True_solution','2003B'); +inputdir.Pivchal03B.truesolbase='solutiononly'; + +%2005B +inputdir.Pivchal05B.pranasol=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_2005B','Case_B','CreateMultiframe','PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.Pivchal05B.MCdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.IMdir=fullfile(basedirname,'Results',runname,'2005B'); +inputdir.Pivchal05B.pranabaseMC='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.pranabaseIM='PIV_scc_chal05_diaCClsq_pass4_'; +inputdir.Pivchal05B.davisbase='B'; +inputdir.Pivchal05B.truesoldir=fullfile(basedirname,'Results','True_solution','2005B'); +inputdir.Pivchal05B.truesolbase='PIVchal05Btruesol'; + +%stagnation_flow +inputdir.stagflow.pranasol=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_stagnationflow','stagnation_images','TR_PIV_MP(1x64x64_75%ov)(3x32x32_50%ov)_nopostprocess'); +inputdir.stagflow.MCdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.IMdir=fullfile(basedirname,'Results',runname,'stagnation_flow'); +inputdir.stagflow.pranabaseMC='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.pranabaseIM='PIV_scc_finalstag_pass4_'; +inputdir.stagflow.davisbase='B'; +inputdir.stagflow.truesoldir=fullfile(basedirname,'Results','True_solution','stagnation_flow'); +inputdir.stagflow.truesolbase='Stagtruesol'; + +%Vortex_Ring +inputdir.vortex.pranasol=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_vortexring','Camera_01','TR_PIV_MP(1x64x64_87%ov)(3x32x32_75%ov)_nopostprocess'); +inputdir.vortex.MCdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.IMdir=fullfile(basedirname,'Results',runname,'Vortex_Ring'); +inputdir.vortex.pranabaseMC='PIV_scc_pass4_'; +inputdir.vortex.pranabaseIM='PIV_scc_pass4_'; +inputdir.vortex.davisbase='B'; +inputdir.vortex.truesoldir=fullfile(basedirname,'Results','True_solution','Vortex_Ring'); +inputdir.vortex.truesolbase='vortextruesol'; + +%Jetdata +inputdir.jet.pranasol=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.davisol=fullfile(basedirname,'uncertainty_tests_Davis','uncertainty_test_jetflow','MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501','TR_PIV_MP(1x32x32_75%ov)(3x16x16_75%ov)_postprocess'); +inputdir.jet.MCdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.IMdir=fullfile(basedirname,'Results',runname,'Jetdata'); +inputdir.jet.pranabaseMC='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.pranabaseIM='PIV_scc_jet_cropped_pass4_'; +inputdir.jet.davisbase='B'; +inputdir.jet.truesoldir=fullfile(basedirname,'Results','True_solution','Jetdata'); +inputdir.jet.truesolbase='jettruesol'; + +end + diff --git a/postprocess_codes/pivchal05b_ixx_postprocess.m b/postprocess_codes/pivchal05b_ixx_postprocess.m new file mode 100644 index 0000000..cfdf865 --- /dev/null +++ b/postprocess_codes/pivchal05b_ixx_postprocess.m @@ -0,0 +1,681 @@ +% uncertainty IXX testing for 2005B case + +clear all; +addpath Z:\Planar_Uncertainty_work\codes\export_fig ; +istring1=sprintf(['%%s%%0%0.0fd.','mat'],3); +istring2=sprintf(['%%s%%0%0.0fd.','txt'],3);%sprintf('%%s%%0%0.0fd',3); +savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\'; + + +% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2005B\'; +% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2005B\newixx\'; +outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\2005B\ws32IM\'; + +biasfact=0; + +if ~exist(outdir1,'dir') + mkdir(outdir1); +end +casename1='cclsq2_'; +outdir=[outdir1,casename1]; +dstore=0; +printfig=0; + +% dirscc='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\2005B\PIV_scc_chal05_diaCClsq_pass4_'; +% dirscc='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\2005B\newtest\PIV_scc_chal05_diaCClsq_pass4_'; +% dirscc2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2005B\PIV_scc_chal05_diaCClsq_pass5_'; + +% dirscc='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\PIV_scc_chal05_diaCClsq_pass4_'; + +% dirscc='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\newixx\PIV_scc_chal05_diaCClsq_pass4_'; +% dirscc2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\PIV_scc_chal05_diaCClsq_pass5_'; +% dirscc='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2005B\PIV_scc_chal05_diaCClsq_pass4_'; +% dirscc2='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2005B\PIV_scc_chal05_diaCClsq_pass5_'; + +% dirscc='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2005B\PIV_scc_chal05_diaCClsq_pass4_'; +% dirscc2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2005B\ws32IM\PIV_scc_chal05_diaCClsq_pass4_'; +dirscc='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2005B\PIV_scc_chal05_diaCClsq_pass4_'; +dirscc2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2005B\PIV_scc_chal05_diaCClsq_pass4_'; + + +% TRUE SOLUTION DIRECTORY +dirtrue='Z:\My_planar_uncertainty_work_and_other_stuff_that_I_need_to_sort\Uncertainty Work\2005B\resultsixx\newresults\CaseBt'; + +f=[9 29 49 69 89 109]; +D=[2.8 2.8]; +toff=0; +for k=1:6 + + n=f(k+toff); + +sccsol=load(sprintf(istring1,dirscc,n)); +sccsol2=load(sprintf(istring1,dirscc2,n)); + +truesol=load(sprintf(istring2,dirtrue,n+1)); + + +%Getting True Solution +x=truesol(:,1); +y=truesol(:,2); +u=truesol(:,3); +v=truesol(:,4); + +Xt=zeros(41,88); +Yt=zeros(41,88); +Ut=zeros(41,88); +Vt=zeros(41,88); +Nx=41;Ny=88; +m=1; +for i=1:41 + for j=1:88 + Xt(i,j)=x(m); + Yt(i,j)=y(m); + Ut(i,j)=u(m); + Vt(i,j)=v(m); + m=m+1; + end +end + +Utrue(:,:,k)=Ut; +Vtrue(:,:,k)=Vt; + +% keyboard; +%Getting Solution + +Us(:,:,k)=(sccsol.U(:,:,1)./2); +Vs(:,:,k)=(sccsol.V(:,:,1)./2); +Eval(:,:,k)=sccsol.Eval(:,:,1); + + +MI(:,:,k)=reshape(sccsol.uncertainty(:,5),Nx,Ny); +IXX2(:,:,k)=reshape(sccsol2.uncertainty(:,9)./2,Nx,Ny); +IYY2(:,:,k)=reshape(sccsol2.uncertainty(:,10)./2,Nx,Ny); +% IXX2(:,:,k)=reshape(sccsol3.uncertainty(:,25)./2,Nx,Ny); +% IYY2(:,:,k)=reshape(sccsol3.uncertainty(:,26)./2,Nx,Ny); + +Autod(:,:,k)=reshape(sccsol.uncertainty(:,6),Nx,Ny); + +Ixx=reshape((sccsol.uncertainty(:,7))./2,Nx,Ny); +Iyy=reshape((sccsol.uncertainty(:,8))./2,Nx,Ny); + +Gpx(:,:,k)=reshape(sccsol.uncertainty(:,15),Nx,Ny); +Gpy(:,:,k)=reshape(sccsol.uncertainty(:,16),Nx,Ny); +Spx(:,:,k)=biasfact*reshape(sccsol.uncertainty(:,17),Nx,Ny); +Spy(:,:,k)=biasfact*reshape(sccsol.uncertainty(:,18),Nx,Ny); + +% %correct gradient (if desired) +[Udiff]=socdiff(Us(:,:,k),16,1); +[Vdiff]=socdiff(Vs(:,:,k),16,2); +% Ixxt = sqrt(Ixx.^2 - (D(1)^2/16)*(Udiff).^2 ); +% Iyyt = sqrt(Iyy.^2 - (D(2)^2/16)*(Vdiff).^2 ); +% Ixxt=Ixx; +% Iyyt=Iyy; +Ixxt= real(sqrt(Ixx.^2 - (Autod(:,:,k).^2/16).*(Udiff).^2 )); +Iyyt= real(sqrt(Iyy.^2 - (Autod(:,:,k).^2/16).*(Vdiff).^2 )); + +%Bias error +mewx(:,:,k)=(abs(Spx(:,:,k)-Gpx(:,:,k)))./2; +mewy(:,:,k)=(abs(Spy(:,:,k)-Gpy(:,:,k)))./2; + +%Scaling for ixx +Ixxt=Ixxt./sqrt(MI(:,:,k)); +Iyyt=Iyyt./sqrt(MI(:,:,k)); +% Final second order moment uncertainty +IXX1(:,:,k)=sqrt(Ixxt.^2+mewx(:,:,k).^2); +IYY1(:,:,k)=sqrt(Iyyt.^2+mewy(:,:,k).^2); +% IXX1(:,:,k)=sqrt(Ixxt.^2); +% IYY1(:,:,k)=sqrt(Iyyt.^2); + +err_sccu(:,:,k)=abs(Utrue(:,:,k)-Us(:,:,k)); +err_sccv(:,:,k)=abs(Vtrue(:,:,k)-Vs(:,:,k)); + +noabserrx(:,:,k)=(Utrue(:,:,k)-Us(:,:,k)); +noabserry(:,:,k)=(Vtrue(:,:,k)-Vs(:,:,k)); + +th=0.2; +tmp_eu=err_sccu(:,:,k); +tmp_ev=err_sccv(:,:,k); +tmp_eu(tmp_eu1); +% % tex(iex)=[]; +% % tixx(iex)=[]; +% +% rms_dx(k)=rms(tex); +% rms_ixx(k)=rms(tixx); + + +% k1=k1+1; +% clear Us Vs Ut Vt MI Ixx Iyy Ixxd Iyyd; +% % figure;imagesc(Us); +% % figure;imagesc(Ut); +% % % figure;imagesc(Vs); +% % % figure;imagesc(Vt); +% % % %figure;hist(MI(:),100); +% keyboard; + +% x_centroid(:,k1)=sccsol.uncertainty(:,21); +% y_centroid(:,k1)=sccsol.uncertainty(:,22); +% x_centroid_scc(:,k1)=sccsol.uncertainty(:,23); +% y_centroid_scc(:,k1)=sccsol.uncertainty(:,24); +if dstore==1 +% Saving Diameters +DCCx(:,:,k)=reshape(sccsol.uncertainty(:,21),Nx,Ny); +DCCy(:,:,k)=reshape(sccsol.uncertainty(:,22),Nx,Ny); +DGccx(:,:,k)=reshape(sccsol.uncertainty(:,29),Nx,Ny); +DGccy(:,:,k)=reshape(sccsol.uncertainty(:,30),Nx,Ny); +Dauto1x(:,:,k)=reshape(sccsol.uncertainty(:,23),Nx,Ny); +Dauto1y(:,:,k)=reshape(sccsol.uncertainty(:,24),Nx,Ny); +Dauto2x(:,:,k)=reshape(sccsol.uncertainty(:,27),Nx,Ny); +Dauto2y(:,:,k)=reshape(sccsol.uncertainty(:,28),Nx,Ny); +D1x(:,:,k)=reshape(sccsol.uncertainty(:,31),Nx,Ny); +D1y(:,:,k)=reshape(sccsol.uncertainty(:,32),Nx,Ny); +end +% keyboard; +end +% keyboard; +% save('2005Brmserror_prana_valid72.mat','rms_dx','rms_dy','rms_ixx','rms_iyy'); +% keyboard; + +% figure;hist(MI(:),100); +% dxz=find(IXX3(:)==-10); +% dyz=find(IYY3(:)==-10); +%[length(dxz) length(dyz)] + +% clim=0.04; +% figure; +% % subplot(2,3,1);imagesc(err_sccu(:,:,1));caxis([0 clim]);title('error at t=9'); +% % subplot(2,3,2);imagesc(err_sccu(:,:,3));caxis([0 clim]);title('error at t=49'); +% % subplot(2,3,3);imagesc(err_sccu(:,:,5));caxis([0 clim]);title('error at t=89'); +% % subplot(2,3,4);imagesc(IXX1(:,:,1));caxis([0 clim]);title('Unx IXX at t=9'); +% % subplot(2,3,5);imagesc(IXX1(:,:,3));caxis([0 clim]);title('Unx IXX at t=49'); +% % subplot(2,3,6);imagesc(IXX1(:,:,5));caxis([0 clim]);title('Unx IXX at t=89'); +% +% subplot(2,3,1);imagesc(err_sccu(:,:,1));caxis([0 clim]);title('error at t=10'); +% subplot(2,3,2);imagesc(err_sccu(:,:,2));caxis([0 clim]);title('error at t=30'); +% subplot(2,3,3);imagesc(err_sccu(:,:,3));caxis([0 clim]);title('error at t=50'); +% subplot(2,3,4);imagesc(err_sccu(:,:,4));caxis([0 clim]);title('error at t=70'); +% subplot(2,3,5);imagesc(err_sccu(:,:,5));caxis([0 clim]);title('error at t=90'); +% subplot(2,3,6);imagesc(err_sccu(:,:,6));caxis([0 clim]);title('error at t=110'); +% keyboard; +% % clim=3; +% % figure; +% % subplot(2,3,1);imagesc(Evalu(:,:,1));caxis([0 clim]);title('error at t=10'); +% % subplot(2,3,2);imagesc(Evalu(:,:,2));caxis([0 clim]);title('error at t=30'); +% % subplot(2,3,3);imagesc(Evalu(:,:,3));caxis([0 clim]);title('error at t=50'); +% % subplot(2,3,4);imagesc(Evalu(:,:,4));caxis([0 clim]);title('error at t=70'); +% % subplot(2,3,5);imagesc(Evalu(:,:,5));caxis([0 clim]);title('error at t=90'); +% % subplot(2,3,6);imagesc(Evalu(:,:,6));caxis([0 clim]);title('error at t=110'); +% % +% % clim=3; +% % figure; +% % subplot(2,3,1);imagesc(Evalv(:,:,1));caxis([0 clim]);title('error at t=10'); +% % subplot(2,3,2);imagesc(Evalv(:,:,2));caxis([0 clim]);title('error at t=30'); +% % subplot(2,3,3);imagesc(Evalv(:,:,3));caxis([0 clim]);title('error at t=50'); +% % subplot(2,3,4);imagesc(Evalv(:,:,4));caxis([0 clim]);title('error at t=70'); +% % subplot(2,3,5);imagesc(Evalv(:,:,5));caxis([0 clim]);title('error at t=90'); +% % subplot(2,3,6);imagesc(Evalv(:,:,6));caxis([0 clim]);title('error at t=110'); + +%% +noabserrx=noabserrx(:); +noabserry=noabserry(:); +err_sccu=err_sccu(:); +err_sccv=err_sccv(:); +IXX1=IXX1(:); +IYY1=IYY1(:); +IXX2=IXX2(:); +IYY2=IYY2(:); + +% IXX2=abs(IXX2); +% IYY2=abs(IYY2); + +MI=MI(:); +if dstore==1 + DCCx=DCCx(:); + DCCy=DCCy(:); + DGccx=DGccx(:); + DGccy=DGccy(:); + Dauto1x=Dauto1x(:); + Dauto1y=Dauto1y(:); + Dauto2x=Dauto2x(:); + Dauto2y=Dauto2y(:); + D1x=D1x(:); + D1y=D1y(:); +end +% DGccx=DGccx./sqrt(MI); +% DGccy=DGccy./sqrt(MI); + +% save([savematdir,'2005BSOM.mat'],'err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2'); + +% save('2005BSOMnew_toterr.mat','err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2'); +% keyboard; +%% +%zero value removal +inx=(find(IXX1(:)==0)); +iny=(find(IYY1(:)==0)); +% keyboard; +IXX1(inx)=[]; +IYY1(iny)=[]; +IXX2(inx)=[]; +IYY2(iny)=[]; +err_sccu(inx)=[]; +err_sccv(iny)=[]; +noabserrx(inx)=[]; +noabserry(iny)=[]; +% MI(inx)=[]; +%% +%nan value removal +inx2=(find(isnan(IXX2))); +iny2=(find(isnan(IXX2))); + +%keyboard; +IXX1(inx2)=[]; +IYY1(iny2)=[]; +IXX2(inx2)=[]; +IYY2(iny2)=[]; +err_sccu(inx2)=[]; +err_sccv(iny2)=[]; + +inx3=(find(isnan(IXX1))); +iny3=(find(isnan(IXX1))); +IXX1(inx3)=[]; +IYY1(iny3)=[]; +IXX2(inx3)=[]; +IYY2(iny3)=[]; +err_sccu(inx3)=[]; +err_sccv(iny3)=[]; + +[length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)] +%% +%invalid vector removal +elm=0.2;%0.2; +inx1=(find(err_sccu(:)>elm)); +iny1=(find(err_sccv(:)>elm)); +[length(inx1) length(iny1)] +%keyboard; +IXX1(inx1)=[]; +IYY1(iny1)=[]; +IXX2(inx1)=[]; +IYY2(iny1)=[]; +err_sccu(inx1)=[]; +err_sccv(iny1)=[]; + +if dstore==1 + DCCx(inx1)=[]; + DGccx(inx1)=[]; + Dauto1x(inx1)=[]; + Dauto2x(inx1)=[]; + D1x(inx1)=[]; + + DCCy(iny1)=[]; + DGccy(iny1)=[]; + Dauto1y(iny1)=[]; + Dauto2y(iny1)=[]; + D1y(iny1)=[]; +end +noabserrx(inx1)=[]; +noabserry(iny1)=[]; +%% +figure;hist(noabserrx,30); +figure;hist(noabserry,30); + +[mean(noabserrx) mean(noabserry)] + +%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%keyboard; + +gmx=(find(err_sccu(:)<=1)); +gmy=(find(err_sccv(:)<=1)); +ngmx=length(gmx); +ngmy=length(gmy); + +%keyboard; +% N=r*s*t; +% [(N-ngmx)*100/N (N-ngmy)*100/N]; + +cnt=0; +for j=1:length(IXX1) +if err_sccu(j)<=IXX1(j) +cnt=cnt+1; +end +end +covx=100*cnt/ngmx; + +cnt=0; +for j=1:length(IYY1) +if err_sccv(j)<=IYY1(j) +cnt=cnt+1; +end +end +covy=100*cnt/ngmy; + +cnt=0; +for j=1:length(IXX2) +if err_sccu(j)<=IXX2(j) +cnt=cnt+1; +end +end +covx2=100*cnt/ngmx; + +cnt=0; +for j=1:length(IYY2) +if err_sccv(j)<=IYY2(j) +cnt=cnt+1; +end +end +covy2=100*cnt/ngmy; + + [covx covx2] + [covy covy2] +% save('2005BSOM_new.mat','err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2'); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% +lw=3; +fs=18; +%% +% % vec=linspace(0,0.2,60); +% % +% % Nrex=histc(err_sccu,vec); +% % Nrey=histc(err_sccv,vec); +% % Nrdx=histc(IXX1,vec); +% % Nrdy=histc(IYY1,vec); +% % Nrdx2=histc(IXX2,vec); +% % Nrdy2=histc(IYY2,vec); +% % +% % ll=1:100:15e2;%max(Nrdx); +% % ll2=1:100:20e2;%max(Nrdy); +% % +% % sigu=rms(err_sccu(err_sccu<=0.15)).*ones(length(ll),1); +% % sigv=rms(err_sccv(err_sccv<=0.15)).*ones(length(ll2),1); +% % +% % mixx1=median(IXX1(IXX1<=0.15)).*ones(length(ll),1); +% % miyy1=median(IYY1(IYY1<=0.15)).*ones(length(ll2),1); +% % +% % mixx2=median(IXX2(IXX2<=0.15)).*ones(length(ll),1); +% % miyy2=median(IYY2(IYY2<=0.15)).*ones(length(ll2),1); +% % +% % +% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); +% % plot(vec,Nrex,'b-o',vec,Nrdx,'r-*',vec,Nrdx2,'m-o',sigu,ll,'k--',mixx1,ll,'r--',mixx2,ll,'m--'); +% % %legend('|errorx|','IXX AMM','IXX DMM','rms|errorx|'); +% % legend('|e_x|','I_X_X','\delta_x','\sigma_{|e_x|}','\sim{I_X_X}','\sim{\delta_x}'); +% % title('Error and Uncertainty histogram'); +% % %title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)])); +% % ylabel('No. of Data Points');xlabel('pixels'); +% % +% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); +% % plot(vec,Nrey,'b-o',vec,Nrdy,'r-*',vec,Nrdy2,'m-o',sigv,ll2,'k--',miyy1,ll2,'r--',miyy2,ll2,'m--'); +% % %legend('|Errory|','IYY AMM','IYY DMM','rms|errory|'); +% % legend('|e_y|','I_Y_Y','\delta_y','\sigma_{|e_y|}','\sim{I_Y_Y}','\sim{\delta_y}'); +% % title('Error and Uncertainty histogram'); +% % %title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)])); +% % ylabel('No. of Data Points');xlabel('pixels'); +%% +thresh=0.1; +vec=linspace(0,thresh,40); +% vec2=linspace(0,0.05,100); +% vec3=linspace(0.03,0.1,100); +% vec4=linspace(0,0.03,100); +% vec5=linspace(0.06,0.1,100); + +Nrex=histc(err_sccu,vec); +Nrey=histc(err_sccv,vec); +Nrdx=histc(IXX1,vec); +Nrdy=histc(IYY1,vec); +Nrdx2=histc(IXX2,vec); +Nrdy2=histc(IYY2,vec); + +ll=1:100:15e2;%max(Nrdx); +ll2=1:100:20e2;%max(Nrdy); +% ll=1:100:15e1;%max(Nrdx); +% ll2=1:100:20e1;%max(Nrdy); + +sigu=rms(err_sccu(err_sccu<=thresh)).*ones(length(ll),1); +sigv=rms(err_sccv(err_sccv<=thresh)).*ones(length(ll2),1); + +% mixx1=median(IXX1(IXX1<=0.1)).*ones(length(ll),1); +% miyy1=median(IYY1(IYY1<=0.1)).*ones(length(ll2),1); +% +% mixx2=median(IXX2(IXX2<=0.1)).*ones(length(ll),1); +% miyy2=median(IYY2(IYY2<=0.1)).*ones(length(ll2),1); + +mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1); +miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + +mixx2=rms(IXX2(IXX2<=thresh)).*ones(length(ll),1); +miyy2=rms(IYY2(IYY2<=thresh)).*ones(length(ll2),1); + +fprintf(['RMSerr X','\n']); +fprintf(num2str([sigu(1) mixx1(1) mixx2(1)],'%02.4f\t')); +fprintf(['\nRMSerr Y','\n']); +fprintf(num2str([sigv(1) miyy1(1) miyy2(1)],'%02.4f\t')); +fprintf('\n'); + + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--'); + +hold off; +%legend('|errorx|','IXX AMM','IXX DMM','rms|errorx|'); +%legend('|e_x|','I_X_X','\delta_x','\sigma_{|e_x|}','\sim{I_X_X}','\sim{\delta_x}'); +% title('Error and Uncertainty histogram(X)'); +ylabel('No. of Data Points');xlabel('pixels'); +box on; +title(sprintf(['Ixx=',num2str(covx),' IMx=',num2str(covx2)])); +if printfig==1 + print(gcf,'-dpng',[outdir,'rmsUhist.png']); +end +% export_fig(gcf,[outdir,'rmsUhist.png']); + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--'); +% plot(vec2,Nppryl,'color',rgb('DarkCyan')); +% plot(vec3,Nppryh,'color',rgb('DarkCyan')); +% plot(vec4,Nmiyl,'color',rgb('SkyBlue')); +% plot(vec5,Nmiyh,'color',rgb('SkyBlue')); +hold off; +%legend('|Errory|','IYY AMM','IYY DMM','rms|errory|'); +%legend('|e_y|','I_Y_Y','\delta_y','\sigma_{|e_y|}','\sim{I_Y_Y}','\sim{\delta_y}'); +% title('Error and Uncertainty histogram(Y)'); +ylabel('No. of Data Points');xlabel('pixels'); +box on; +title(sprintf(['Iyy=',num2str(covy),' IMy=',num2str(covy2)])); +if printfig==1 + print(gcf,'-dpng',[outdir,'rmsVhist.png']); +end + %% +Nbin=20; + +maxixx=0.08; +minixx=0; +vu=linspace(minixx,maxixx,Nbin); +maxiyy=0.08; +miniyy=0; +vv=linspace(miniyy,maxiyy,Nbin); + +rmserru1=zeros(1,length(vu)); +rmserrv1=zeros(1,length(vv)); +getixx1=zeros(1,length(vu)); +getiyy1=zeros(1,length(vv)); +numcx1=zeros(1,length(vu)); +numcy1=zeros(1,length(vv)); +covxbin1=zeros(1,length(vu)); +covybin1=zeros(1,length(vv)); + +rmserru2=zeros(1,length(vu)); +rmserrv2=zeros(1,length(vv)); +getixx2=zeros(1,length(vu)); +getiyy2=zeros(1,length(vv)); +numcx2=zeros(1,length(vu)); +numcy2=zeros(1,length(vv)); +covxbin2=zeros(1,length(vu)); +covybin2=zeros(1,length(vv)); + +%p=1; + +for p=1:length(vu) + if p=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); + + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); + + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1veccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; + %% +% err_up=err_up-mean(err_up); +% err_vp=err_vp-mean(err_vp); + %% Calculate coverage + ngmx1=length(err_up); + ngmy1=length(err_vp); + ngmx2=length(err_ud); + ngmy2=length(err_vd); + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UMCx(j) && err_up(j)>-UMCx(j) + cnt=cnt+1; + end + end + covxMC=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UMCy(j) && err_vp(j)>-UMCy(j) + cnt=cnt+1; + end + end + covyMC=100*cnt/ngmy1; + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UIMx(j) && err_up(j)>-UIMx(j) + cnt=cnt+1; + end + end + covxIM=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UIMy(j) && err_vp(j)>-UIMy(j) + cnt=cnt+1; + end + end + covyIM=100*cnt/ngmy1; + + + cnt=0; + for j=1:ngmx2 + if err_ud(j)<=UCSx(j) && err_ud(j)>-UCSx(j) + cnt=cnt+1; + end + end + covxCS=100*cnt/ngmx2; + + cnt=0; + for j=1:ngmy2 + if err_vd(j)<=UCSy(j) && err_vd(j)>-UCSy(j) + cnt=cnt+1; + end + end + covyCS=100*cnt/ngmy2; + + + [covxMC covyMC covxIM covyIM covxCS covyCS] +% [rms(err_up) rms(err_vp) rms(err_ud) rms(err_vd)] +% keyboard; + %% Plot Functions call + %Plot X and Y error and uncertainty histogram together for all 3 metrics + + error_prana=[err_up;err_vp]; + error_davis=[err_ud;err_vd]; + UMC=[UMCx;UMCy]; + UIM=[UIMx;UIMy]; + UCS=[UCSx;UCSy]; + %% Find percentage of error within Std deviation error + cutoff=0.1;%0.1%0.5 + %Prana + eP=error_prana; + eP(eP>cutoff)=[]; + eP=eP-mean(eP); + ePstd=std(eP); + indexP=find(abs(eP)<=ePstd); + covP=100*length(indexP)/length(eP); + %DaVis + eD=error_davis; + eD(eD>cutoff)=[]; + eD=eD-mean(eD); + eDstd=std(eD); + indexD=find(abs(eD)<=eDstd); + covD=100*length(indexD)/length(eD); + + [covP covD] + + %% + Nbin1=40; + ler=-0.1;uer=0.1; + Nbin2=60; + lun=0;uun=0.1; + + vec=linspace(ler,uer,Nbin1); + Nep=histc(error_prana,vec); + Ned=histc(error_davis,vec); + + vec2=linspace(lun,uun,Nbin2); + Nmc=histc(UMC,vec2); + Nim=histc(UIM,vec2); + Ncs=histc(UCS,vec2); + + Lp=length(error_prana); + Ld=length(error_davis); + Nep=(Nep/Lp)*100; + Ned=(Ned/Ld)*100; + Nmc=(Nmc/Lp)*100; + Nim=(Nim/Lp)*100; + Ncs=(Ncs/Ld)*100; + + + lw=1.2;fs=20; + set(0,'DefaultAxesFontName', 'Times New Roman'); + crimson=rgb('crimson'); + dblue=rgb('deepskyblue'); + blueviolet=rgb('blueviolet'); + figure; + set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(vec,Nep,'k-'); + plot(vec,Ned,'k--'); + plot(vec2,Nmc,'r-','MarkerSize',4,'color',crimson); + plot(vec2,Nim,'c-','MarkerSize',4,'color',dblue); + plot(vec2,Ncs,'g-','MarkerSize',4,'color',blueviolet); + + yl=ylim(gca); + + ll=linspace(0,maxyl,10); +% ll2=1:100:max(0.75*yl(2)); + sigprana=rms(error_prana(abs(error_prana)<=0.1)).*ones(length(ll),1); + sigdavis=rms(error_davis(abs(error_davis)<=0.1)).*ones(length(ll),1); + sigMC=rms(UMC(UMC<=0.1)).*ones(length(ll),1); + sigIM=rms(UIM(UIM<=0.1)).*ones(length(ll),1); + sigCS=rms(UCS(UCS<=0.1)).*ones(length(ll),1); + + [mean(error_prana) mean(error_davis)] + + plot(sigprana,ll,'k-'); + plot(sigdavis,ll,'k--'); + plot(sigMC,ll,'r-','MarkerSize',4,'color',crimson); + plot(sigIM,ll,'c-.','MarkerSize',4,'color',dblue); + plot(sigCS,ll,'g--','MarkerSize',4,'color',blueviolet); + + hold off; +% ylabel('No. of Data Points'); +% xlabel('Error and Uncertainty (pixels)'); +% lbox=legend('e_{prana}','e_{davis}','MC','PD','CS','location','northeast'); +% set(lbox,'FontSize',16); +% title('Turbulent Boundary Layer \newline'); + grid on;box on; + set(gcf,'color','white'); +% axis square; + hold off; + if printfig==1 + export_fig(gcf,fullfile(outputdir{2},[casename{1},'_Histogram.png']),'-painters','-r360'); + end + if savematplot==1 + save(fullfile(outputdir{1},'PivChal03B_hist.mat'),'vec','Nep','Ned','vec2','Nmc','Nim','Ncs','ll','sigprana','sigdavis','sigMC','sigIM','sigCS'); + end + + [sigprana(1) sigdavis(1) sigMC(1) sigIM(1) sigCS(1)] + + %% Plot RMS error in uncertainty bins scatter plot + % For MC + Nbin3=8; + Ubin1=linspace(0, 0.08, Nbin3); + + numcx1=zeros(Nbin3,1); + covxbin1=zeros(Nbin3,1); + rmserru1=zeros(Nbin3,1); + rmsMC=zeros(Nbin3,1); + + for p=1:length(Ubin1) + if p=Ubin1(p))); + elseif p==length(Ubin1) + [val1]= find((UMC>=Ubin1(p))); + end + numcx1(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UMC(val1); + cvrx1=find(tmperr=Ubin2(p))); + elseif p==length(Ubin2) + [val1]= find((UIM>=Ubin2(p))); + end + numcx2(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UIM(val1); + cvrx1=find(tmperr=Ubin3(p))); + elseif p==length(Ubin3) + [val1]= find((UCS>=Ubin3(p))); + end + numcx3(p)=(length(val1)); + tmperr=(error_davis(val1)); + tempMC=UCS(val1); + cvrx1=find(tmperrveccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; + + %% Calculate coverage + ngmx1=length(err_up); + ngmy1=length(err_vp); + ngmx2=length(err_ud); + ngmy2=length(err_vd); + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UMCx(j) && err_up(j)>-UMCx(j) + cnt=cnt+1; + end + end + covxMC=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UMCy(j) && err_vp(j)>-UMCy(j) + cnt=cnt+1; + end + end + covyMC=100*cnt/ngmy1; + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UIMx(j) && err_up(j)>-UIMx(j) + cnt=cnt+1; + end + end + covxIM=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UIMy(j) && err_vp(j)>-UIMy(j) + cnt=cnt+1; + end + end + covyIM=100*cnt/ngmy1; + + + cnt=0; + for j=1:ngmx2 + if err_ud(j)<=UCSx(j) && err_ud(j)>-UCSx(j) + cnt=cnt+1; + end + end + covxCS=100*cnt/ngmx2; + + cnt=0; + for j=1:ngmy2 + if err_vd(j)<=UCSy(j) && err_vd(j)>-UCSy(j) + cnt=cnt+1; + end + end + covyCS=100*cnt/ngmy2; + + + [covxMC covyMC covxIM covyIM covxCS covyCS] +% [rms(err_up) rms(err_vp) rms(err_ud) rms(err_vd)] +% keyboard; + %% Plot Functions call + %Plot X and Y error and uncertainty histogram together for all 3 metrics + + error_prana=[err_up;err_vp]; + error_davis=[err_ud;err_vd]; + UMC=[UMCx;UMCy]; + UIM=[UIMx;UIMy]; + UCS=[UCSx;UCSy]; + + %% Find percentage of error within Std deviation error + cutoff=0.1;%0.1 + %Prana + eP=error_prana; + eP(eP>cutoff)=[]; + eP=eP-mean(eP); + ePstd=std(eP); + indexP=find(abs(eP)<=ePstd); + covP=100*length(indexP)/length(eP); + %DaVis + eD=error_davis; + eD(eD>cutoff)=[]; + eD=eD-mean(eD); + eDstd=std(eD); + indexD=find(abs(eD)<=eDstd); + covD=100*length(indexD)/length(eD); + + [covP covD] + + %% + Nbin1=40; + ler=-0.1;uer=0.1; + Nbin2=60; + lun=0;uun=0.1; + + vec=linspace(ler,uer,Nbin1); + Nep=histc(error_prana,vec); + Ned=histc(error_davis,vec); + + vec2=linspace(lun,uun,Nbin2); + Nmc=histc(UMC,vec2); + Nim=histc(UIM,vec2); + Ncs=histc(UCS,vec2); + + Lp=length(error_prana); + Ld=length(error_davis); + Nep=(Nep/Lp)*100; + Ned=(Ned/Ld)*100; + Nmc=(Nmc/Lp)*100; + Nim=(Nim/Lp)*100; + Ncs=(Ncs/Ld)*100; + + lw=1.2;fs=20; + set(0,'DefaultAxesFontName', 'Times New Roman'); + crimson=rgb('crimson'); + dblue=rgb('deepskyblue'); + blueviolet=rgb('blueviolet'); + figure; + set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(vec,Nep,'k-'); + plot(vec,Ned,'k--'); + plot(vec2,Nmc,'r-','MarkerSize',4,'color',crimson); + plot(vec2,Nim,'c-','MarkerSize',4,'color',dblue); + plot(vec2,Ncs,'g-','MarkerSize',4,'color',blueviolet); + + yl=ylim(gca); + + ll=linspace(0,maxyl,10); +% ll2=1:100:max(0.75*yl(2)); + sigprana=rms(error_prana(abs(error_prana)<=0.1)).*ones(length(ll),1); + sigdavis=rms(error_davis(abs(error_davis)<=0.1)).*ones(length(ll),1); + sigMC=rms(UMC(UMC<=0.1)).*ones(length(ll),1); + sigIM=rms(UIM(UIM<=0.1)).*ones(length(ll),1); + sigCS=rms(UCS(UCS<=0.1)).*ones(length(ll),1); + + [mean(error_prana) mean(error_davis)] + + plot(sigprana,ll,'k-'); + plot(sigdavis,ll,'k--'); + plot(sigMC,ll,'r-','MarkerSize',4,'color',crimson); + plot(sigIM,ll,'c-.','MarkerSize',4,'color',dblue); + plot(sigCS,ll,'g--','MarkerSize',4,'color',blueviolet); + + hold off; +% ylabel('No. of Data Points'); +% xlabel('Error and Uncertainty (pixels)'); +% lbox=legend('e_{prana}','e_{davis}','MC','PD','CS','location','northeast'); +% set(lbox,'FontSize',16); +% title('Laminar Separation Bubble \newline'); + grid on;box on; + set(gcf,'color','white'); +% axis square; + hold off; + if printfig==1 + export_fig(gcf,fullfile(outputdir{2},[casename{2},'_Histogram.png']),'-painters','-r360'); + end + if savematplot==1 + save(fullfile(outputdir{1},'PivChal05B_hist.mat'),'vec','Nep','Ned','vec2','Nmc','Nim','Ncs','ll','sigprana','sigdavis','sigMC','sigIM','sigCS'); + end + [sigprana(1) sigdavis(1) sigMC(1) sigIM(1) sigCS(1)] + + %% Plot RMS error in uncertainty bins scatter plot + + % For MC + Nbin3=8; + Ubin1=linspace(0, 0.08, Nbin3); + + numcx1=zeros(Nbin3,1); + covxbin1=zeros(Nbin3,1); + rmserru1=zeros(Nbin3,1); + rmsMC=zeros(Nbin3,1); + + for p=1:length(Ubin1) + if p=Ubin1(p))); + elseif p==length(Ubin1) + [val1]= find((UMC>=Ubin1(p))); + end + numcx1(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UMC(val1); + cvrx1=find(tmperr=Ubin2(p))); + elseif p==length(Ubin2) + [val1]= find((UIM>=Ubin2(p))); + end + numcx2(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UIM(val1); + cvrx1=find(tmperr=Ubin3(p))); + elseif p==length(Ubin3) + [val1]= find((UCS>=Ubin3(p))); + end + numcx3(p)=(length(val1)); + tmperr=(error_davis(val1)); + tempMC=UCS(val1); + cvrx1=find(tmperrveccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; + + %% Calculate coverage + ngmx1=length(err_up); + ngmy1=length(err_vp); + ngmx2=length(err_ud); + ngmy2=length(err_vd); + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UMCx(j) && err_up(j)>-UMCx(j) + cnt=cnt+1; + end + end + covxMC=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UMCy(j) && err_vp(j)>-UMCy(j) + cnt=cnt+1; + end + end + covyMC=100*cnt/ngmy1; + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UIMx(j) && err_up(j)>-UIMx(j) + cnt=cnt+1; + end + end + covxIM=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UIMy(j) && err_vp(j)>-UIMy(j) + cnt=cnt+1; + end + end + covyIM=100*cnt/ngmy1; + + + cnt=0; + for j=1:ngmx2 + if err_ud(j)<=UCSx(j) && err_ud(j)>-UCSx(j) + cnt=cnt+1; + end + end + covxCS=100*cnt/ngmx2; + + cnt=0; + for j=1:ngmy2 + if err_vd(j)<=UCSy(j) && err_vd(j)>-UCSy(j) + cnt=cnt+1; + end + end + covyCS=100*cnt/ngmy2; + + + [covxMC covyMC covxIM covyIM covxCS covyCS] +% [rms(err_up) rms(err_vp) rms(err_ud) rms(err_vd)] +% keyboard; + %% Plot Functions call + %Plot X and Y error and uncertainty histogram together for all 3 metrics + + error_prana=[err_up;err_vp]; + error_davis=[err_ud;err_vd]; + UMC=[UMCx;UMCy]; + UIM=[UIMx;UIMy]; + UCS=[UCSx;UCSy]; + + %% Find percentage of error within Std deviation error + cutoff=0.2;%0.2 + %Prana + eP=error_prana; + eP(eP>cutoff)=[]; + eP=eP-mean(eP); + ePstd=std(eP); + indexP=find(abs(eP)<=ePstd); + covP=100*length(indexP)/length(eP); + %DaVis + eD=error_davis; + eD(eD>cutoff)=[]; + eD=eD-mean(eD); + eDstd=std(eD); + indexD=find(abs(eD)<=eDstd); + covD=100*length(indexD)/length(eD); + + [covP covD] + + %% + Nbin1=40; + ler=-0.3;uer=0.3; + Nbin2=60; + lun=0;uun=0.3; + + vec=linspace(ler,uer,Nbin1); + Nep=histc(error_prana,vec); + Ned=histc(error_davis,vec); + + vec2=linspace(lun,uun,Nbin2); + Nmc=histc(UMC,vec2); + Nim=histc(UIM,vec2); + Ncs=histc(UCS,vec2); + + Lp=length(error_prana); + Ld=length(error_davis); + Nep=(Nep/Lp)*100; + Ned=(Ned/Ld)*100; + Nmc=(Nmc/Lp)*100; + Nim=(Nim/Lp)*100; + Ncs=(Ncs/Ld)*100; + + lw=1.2;fs=20; + set(0,'DefaultAxesFontName', 'Times New Roman'); + crimson=rgb('crimson'); + dblue=rgb('deepskyblue'); + blueviolet=rgb('blueviolet'); + figure; + set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(vec,Nep,'k-'); + plot(vec,Ned,'k--'); + plot(vec2,Nmc,'r-','MarkerSize',4,'color',crimson); + plot(vec2,Nim,'c-','MarkerSize',4,'color',dblue); + plot(vec2,Ncs,'g-','MarkerSize',4,'color',blueviolet); + + yl=ylim(gca); + + ll=linspace(0,maxyl,10); +% ll2=1:100:max(0.75*yl(2)); + sigprana=rms(error_prana(abs(error_prana)<=0.2)).*ones(length(ll),1); + sigdavis=rms(error_davis(abs(error_davis)<=0.2)).*ones(length(ll),1); + sigMC=rms(UMC(UMC<=0.2)).*ones(length(ll),1); + sigIM=rms(UIM(UIM<=0.2)).*ones(length(ll),1); + sigCS=rms(UCS(UCS<=0.2)).*ones(length(ll),1); + + [mean(error_prana) mean(error_davis)] + + plot(sigprana,ll,'k-'); + plot(sigdavis,ll,'k--'); + plot(sigMC,ll,'r-','MarkerSize',4,'color',crimson); + plot(sigIM,ll,'c-.','MarkerSize',4,'color',dblue); + plot(sigCS,ll,'g--','MarkerSize',4,'color',blueviolet); + + hold off; +% ylabel('No. of Data Points'); +% xlabel('Error and Uncertainty (pixels)'); +% lbox=legend('e_{prana}','e_{davis}','MC','PD','CS','location','northeast'); +% set(lbox,'FontSize',16); +% title('Stagnation Flow \newline'); + grid on;box on; + set(gcf,'color','white'); +% axis square; + hold off; + if printfig==1 + export_fig(gcf,fullfile(outputdir{2},[casename{3},'_Histogram.png']),'-painters','-r360'); + end + if savematplot==1 + save(fullfile(outputdir{1},'stagnation_flow_hist.mat'),'vec','Nep','Ned','vec2','Nmc','Nim','Ncs','ll','sigprana','sigdavis','sigMC','sigIM','sigCS'); + end + [sigprana(1) sigdavis(1) sigMC(1) sigIM(1) sigCS(1)] + + %% Plot RMS error in uncertainty bins scatter plot + + % For MC + Nbin3=8; + Ubin1=linspace(0, 0.3, Nbin3); + + numcx1=zeros(Nbin3,1); + covxbin1=zeros(Nbin3,1); + rmserru1=zeros(Nbin3,1); + rmsMC=zeros(Nbin3,1); + + for p=1:length(Ubin1) + if p=Ubin1(p))); + elseif p==length(Ubin1) + [val1]= find((UMC>=Ubin1(p))); + end + numcx1(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UMC(val1); + cvrx1=find(tmperr=Ubin2(p))); + elseif p==length(Ubin2) + [val1]= find((UIM>=Ubin2(p))); + end + numcx2(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UIM(val1); + cvrx1=find(tmperr=Ubin3(p))); + elseif p==length(Ubin3) + [val1]= find((UCS>=Ubin3(p))); + end + numcx3(p)=(length(val1)); + tmperr=(error_davis(val1)); + tempMC=UCS(val1); + cvrx1=find(tmperrveccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; + + %% Calculate coverage + ngmx1=length(err_up); + ngmy1=length(err_vp); + ngmx2=length(err_ud); + ngmy2=length(err_vd); + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UMCx(j) && err_up(j)>-UMCx(j) + cnt=cnt+1; + end + end + covxMC=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UMCy(j) && err_vp(j)>-UMCy(j) + cnt=cnt+1; + end + end + covyMC=100*cnt/ngmy1; + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UIMx(j) && err_up(j)>-UIMx(j) + cnt=cnt+1; + end + end + covxIM=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UIMy(j) && err_vp(j)>-UIMy(j) + cnt=cnt+1; + end + end + covyIM=100*cnt/ngmy1; + + + cnt=0; + for j=1:ngmx2 + if err_ud(j)<=UCSx(j) && err_ud(j)>-UCSx(j) + cnt=cnt+1; + end + end + covxCS=100*cnt/ngmx2; + + cnt=0; + for j=1:ngmy2 + if err_vd(j)<=UCSy(j) && err_vd(j)>-UCSy(j) + cnt=cnt+1; + end + end + covyCS=100*cnt/ngmy2; + + + [covxMC covyMC covxIM covyIM covxCS covyCS] +% [rms(err_up) rms(err_vp) rms(err_ud) rms(err_vd)] +% keyboard; + %% Plot Functions call + %Plot X and Y error and uncertainty histogram together for all 3 metrics + + error_prana=[err_up;err_vp]; + error_davis=[err_ud;err_vd]; + UMC=[UMCx;UMCy]; + UIM=[UIMx;UIMy]; + UCS=[UCSx;UCSy]; + + %% Find percentage of error within Std deviation error + cutoff=0.2;%0.2 + %Prana + eP=error_prana; + eP(eP>cutoff)=[]; + eP=eP-mean(eP); + ePstd=std(eP); + indexP=find(abs(eP)<=ePstd); + covP=100*length(indexP)/length(eP); + %DaVis + eD=error_davis; + eD(eD>cutoff)=[]; + eD=eD-mean(eD); + eDstd=std(eD); + indexD=find(abs(eD)<=eDstd); + covD=100*length(indexD)/length(eD); + + [covP covD] + + %% + Nbin1=40; + ler=-0.2;uer=0.2; + Nbin2=60; + lun=0;uun=0.2; + + vec=linspace(ler,uer,Nbin1); + Nep=histc(error_prana,vec); + Ned=histc(error_davis,vec); + + vec2=linspace(lun,uun,Nbin2); + Nmc=histc(UMC,vec2); + Nim=histc(UIM,vec2); + Ncs=histc(UCS,vec2); + + Lp=length(error_prana); + Ld=length(error_davis); + Nep=(Nep/Lp)*100; + Ned=(Ned/Ld)*100; + Nmc=(Nmc/Lp)*100; + Nim=(Nim/Lp)*100; + Ncs=(Ncs/Ld)*100; + + lw=1.2;fs=20; + set(0,'DefaultAxesFontName', 'Times New Roman'); + crimson=rgb('crimson'); + dblue=rgb('deepskyblue'); + blueviolet=rgb('blueviolet'); + figure; + set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(vec,Nep,'k-'); + plot(vec,Ned,'k--'); + plot(vec2,Nmc,'r-','MarkerSize',4,'color',crimson); + plot(vec2,Nim,'c-','MarkerSize',4,'color',dblue); + plot(vec2,Ncs,'g-','MarkerSize',4,'color',blueviolet); + + yl=ylim(gca); + + ll=linspace(0,maxyl,10); +% ll2=1:100:max(0.75*yl(2)); + sigprana=rms(error_prana(abs(error_prana)<=0.2)).*ones(length(ll),1); + sigdavis=rms(error_davis(abs(error_davis)<=0.2)).*ones(length(ll),1); + sigMC=rms(UMC(UMC<=0.2)).*ones(length(ll),1); + sigIM=rms(UIM(UIM<=0.2)).*ones(length(ll),1); + sigCS=rms(UCS(UCS<=0.2)).*ones(length(ll),1); + + [mean(error_prana(abs(error_prana)<=0.2)) mean(error_davis(abs(error_davis)<=0.2))] + + plot(sigprana,ll,'k-'); + plot(sigdavis,ll,'k--'); + plot(sigMC,ll,'r-','MarkerSize',4,'color',crimson); + plot(sigIM,ll,'c-.','MarkerSize',4,'color',dblue); + plot(sigCS,ll,'g--','MarkerSize',4,'color',blueviolet); + + hold off; +% ylabel('No. of Data Points'); +% xlabel('Error and Uncertainty (pixels)'); +% lbox=legend('e_{prana}','e_{davis}','MC','PD','CS','location','northeast'); +% set(lbox,'FontSize',16); +% title('Vortex Ring \newline'); + grid on;box on; + set(gcf,'color','white'); +% axis square; + hold off; + if printfig==1 + export_fig(gcf,fullfile(outputdir{2},[casename{4},'_Histogram.png']),'-painters','-r360'); + end + if savematplot==1 + save(fullfile(outputdir{1},'Vortex_Ring_hist.mat'),'vec','Nep','Ned','vec2','Nmc','Nim','Ncs','ll','sigprana','sigdavis','sigMC','sigIM','sigCS'); + end + [sigprana(1) sigdavis(1) sigMC(1) sigIM(1) sigCS(1)] + + %% Plot RMS error in uncertainty bins scatter plot + + % For MC + Nbin3=8; + Ubin1=linspace(0, 0.1, Nbin3); + + numcx1=zeros(Nbin3,1); + covxbin1=zeros(Nbin3,1); + rmserru1=zeros(Nbin3,1); + rmsMC=zeros(Nbin3,1); + + for p=1:length(Ubin1) + if p=Ubin1(p))); + elseif p==length(Ubin1) + [val1]= find((UMC>=Ubin1(p))); + end + numcx1(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UMC(val1); + cvrx1=find(tmperr=Ubin2(p))); + elseif p==length(Ubin2) + [val1]= find((UIM>=Ubin2(p))); + end + numcx2(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UIM(val1); + cvrx1=find(tmperr=Ubin3(p))); + elseif p==length(Ubin3) + [val1]= find((UCS>=Ubin3(p))); + end + numcx3(p)=(length(val1)); + tmperr=(error_davis(val1)); + tempMC=UCS(val1); + cvrx1=find(tmperrveccutoff)); + iny1=(find(abs(err_vp)>veccutoff)); + + inx2=(find(abs(err_ud)>veccutoff)); + iny2=(find(abs(err_vd)>veccutoff)); + + err_up(inx1)=[]; + err_vp(iny1)=[]; + UMCx(inx1)=[]; + UMCy(iny1)=[]; + UIMx(inx1)=[]; + UIMy(iny1)=[]; + + err_ud(inx2)=[]; + err_vd(iny2)=[]; + UCSx(inx2)=[]; + UCSy(iny2)=[]; + + %% Calculate coverage + ngmx1=length(err_up); + ngmy1=length(err_vp); + ngmx2=length(err_ud); + ngmy2=length(err_vd); + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UMCx(j) && err_up(j)>-UMCx(j) + cnt=cnt+1; + end + end + covxMC=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UMCy(j) && err_vp(j)>-UMCy(j) + cnt=cnt+1; + end + end + covyMC=100*cnt/ngmy1; + + cnt=0; + for j=1:ngmx1 + if err_up(j)<=UIMx(j) && err_up(j)>-UIMx(j) + cnt=cnt+1; + end + end + covxIM=100*cnt/ngmx1; + + cnt=0; + for j=1:ngmy1 + if err_vp(j)<=UIMy(j) && err_vp(j)>-UIMy(j) + cnt=cnt+1; + end + end + covyIM=100*cnt/ngmy1; + + + cnt=0; + for j=1:ngmx2 + if err_ud(j)<=UCSx(j) && err_ud(j)>-UCSx(j) + cnt=cnt+1; + end + end + covxCS=100*cnt/ngmx2; + + cnt=0; + for j=1:ngmy2 + if err_vd(j)<=UCSy(j) && err_vd(j)>-UCSy(j) + cnt=cnt+1; + end + end + covyCS=100*cnt/ngmy2; + + + [covxMC covyMC covxIM covyIM covxCS covyCS] +% [rms(err_up) rms(err_vp) rms(err_ud) rms(err_vd)] +% keyboard; + %% Plot Functions call + %Plot X and Y error and uncertainty histogram together for all 3 metrics + + error_prana=[err_up;err_vp]; + error_davis=[err_ud;err_vd]; + UMC=[UMCx;UMCy]; + UIM=[UIMx;UIMy]; + UCS=[UCSx;UCSy]; + + %% Find percentage of error within Std deviation error + cutoff=0.2;%0.2 + %Prana + eP=error_prana; + eP(eP>cutoff)=[]; + eP=eP-mean(eP); + ePstd=std(eP); + indexP=find(abs(eP)<=ePstd); + covP=100*length(indexP)/length(eP); + %DaVis + eD=error_davis; + eD(eD>cutoff)=[]; + eD=eD-mean(eD); + eDstd=std(eD); + indexD=find(abs(eD)<=eDstd); + covD=100*length(indexD)/length(eD); + + [covP covD] + + %% + Nbin1=40; + ler=-0.4;uer=0.4; + Nbin2=60; + lun=0;uun=0.4; + + vec=linspace(ler,uer,Nbin1); + Nep=histc(error_prana,vec); + Ned=histc(error_davis,vec); + + vec2=linspace(lun,uun,Nbin2); + Nmc=histc(UMC,vec2); + Nim=histc(UIM,vec2); + Ncs=histc(UCS,vec2); + + Lp=length(error_prana); + Ld=length(error_davis); + Nep=(Nep/Lp)*100; + Ned=(Ned/Ld)*100; + Nmc=(Nmc/Lp)*100; + Nim=(Nim/Lp)*100; + Ncs=(Ncs/Ld)*100; + + lw=1.2;fs=20; + set(0,'DefaultAxesFontName', 'Times New Roman'); + crimson=rgb('crimson'); + dblue=rgb('deepskyblue'); + blueviolet=rgb('blueviolet'); + figure; + set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); + hold on; + plot(vec,Nep,'k-'); + plot(vec,Ned,'k--'); + plot(vec2,Nmc,'r-','MarkerSize',4,'color',crimson); + plot(vec2,Nim,'c-','MarkerSize',4,'color',dblue); + plot(vec2,Ncs,'g-','MarkerSize',4,'color',blueviolet); + + yl=ylim(gca); + + ll=linspace(0,maxyl,10); +% ll2=1:100:max(0.75*yl(2)); + sigprana=rms(error_prana(abs(error_prana)<=0.2)).*ones(length(ll),1); + sigdavis=rms(error_davis(abs(error_davis)<=0.2)).*ones(length(ll),1); + sigMC=rms(UMC(UMC<=0.2)).*ones(length(ll),1); + sigIM=rms(UIM(UIM<=0.2)).*ones(length(ll),1); + sigCS=rms(UCS(UCS<=0.2)).*ones(length(ll),1); + + [mean(error_prana(abs(error_prana)<=0.2)) mean(error_davis(abs(error_davis)<=0.2))] + + plot(sigprana,ll,'k-'); + plot(sigdavis,ll,'k--'); + plot(sigMC,ll,'r-','MarkerSize',4,'color',crimson); + plot(sigIM,ll,'c-.','MarkerSize',4,'color',dblue); + plot(sigCS,ll,'g--','MarkerSize',4,'color',blueviolet); + + hold off; +% ylabel('No. of Data Points'); +% xlabel('Error and Uncertainty (pixels)'); +% lbox=legend('e_{prana}','e_{davis}','MC','PD','CS','location','northeast'); +% set(lbox,'FontSize',16); +% title('Jet Inviscid Core \newline'); + grid on;box on; + set(gcf,'color','white'); +% axis square; + hold off; + if printfig==1 + export_fig(gcf,fullfile(outputdir{2},[casename{5},'_Histogram.png']),'-painters','-r360'); + end + if savematplot==1 + save(fullfile(outputdir{1},'Jetdata_hist.mat'),'vec','Nep','Ned','vec2','Nmc','Nim','Ncs','ll','sigprana','sigdavis','sigMC','sigIM','sigCS'); + end + [sigprana(1) sigdavis(1) sigMC(1) sigIM(1) sigCS(1)] + + %% Plot RMS error in uncertainty bins scatter plot + + % For MC + Nbin3=8; + Ubin1=linspace(0, 0.12, Nbin3); + + numcx1=zeros(Nbin3,1); + covxbin1=zeros(Nbin3,1); + rmserru1=zeros(Nbin3,1); + rmsMC=zeros(Nbin3,1); + + for p=1:length(Ubin1) + if p=Ubin1(p))); + elseif p==length(Ubin1) + [val1]= find((UMC>=Ubin1(p))); + end + numcx1(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UMC(val1); + cvrx1=find(tmperr=Ubin2(p))); + elseif p==length(Ubin2) + [val1]= find((UIM>=Ubin2(p))); + end + numcx2(p)=(length(val1)); + tmperr=(error_prana(val1)); + tempMC=UIM(val1); + cvrx1=find(tmperr=Ubin3(p))); + elseif p==length(Ubin3) + [val1]= find((UCS>=Ubin3(p))); + end + numcx3(p)=(length(val1)); + tmperr=(error_davis(val1)); + tempMC=UCS(val1); + cvrx1=find(tmperrthresh)); +iny1=(find(err_sccv(:)>thresh)); +%keyboard; +IXX1(inx1)=[]; +IYY1(iny1)=[]; +IXX2(inx1)=[]; +IYY2(iny1)=[]; +err_sccu(inx1)=[]; +err_sccv(iny1)=[]; +[length(inx1) length(iny1)] +% upprx(inx1,:)=[];uprmsx(inx1,:)=[];upcex(inx1,:)=[];uenx(inx1,:)=[];umix(inx1,:)=[]; +% uppry(iny1,:)=[];uprmsy(iny1,:)=[];upcey(iny1,:)=[];ueny(iny1,:)=[];umiy(iny1,:)=[]; + +if dstore==1 + DCCx(inx1)=[]; + DGccx(inx1)=[]; + Dauto1x(inx1)=[]; + Dauto2x(inx1)=[]; + D1x(inx1)=[]; + + DCCy(iny1)=[]; + DGccy(iny1)=[]; + Dauto1y(iny1)=[]; + Dauto2y(iny1)=[]; + D1y(iny1)=[]; +end +noabserrx(inx1)=[]; +noabserry(iny1)=[]; +%% +ll1=-0.4;lu1=0.4; +ll2=-0.2;lu2=0.2; +inx5=find(noabserrx>lu1 | noabserrxlu2 | noabserry=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); +% tmp1(tmp1>1)=[]; +% tmp2(tmp2>1)=[]; + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); +% tmp1(tmp1>1)=[]; +% tmp2(tmp2>1)=[]; + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1=0.01); +% % lny=find(abs(diffy)>=0.01); +% % +% % fprintf([casename,' \n ']); +% % fprintf('Difference in UMCx UMCy between SOCdiff and Richdiff: \n %d %d \n',[length(lnx),length(lny)]); +% % fprintf('Percent Difference in UMCx UMCy between SOCdiff and Richdiff: \n %f %f \n',[100*length(lnx)/length(diffx) 100*length(lny)/length(diffx)]); +% % [mean(diffx(lnx)) mean(diffy(lny))] +% % length(diffx) +% % +% % Udiff1=A1.Udiff; +% % Vdiff1=A1.Vdiff; +% % L11=[length(find(abs(Udiff1)>0.05)) length(find(abs(Vdiff1)>0.05))]; +% % Udiff2=B1.Udiff; +% % Vdiff2=B1.Vdiff; +% % L22=[length(find(abs(Udiff2)>0.05)) length(find(abs(Vdiff2)>0.05))]; +% % +% % fprintf('Udiff and Vdiff >0.05 SOCdiff: \n %d %d \n',L11); +% % fprintf('Udiff and Vdiff >0.05 Richdiff: \n %d %d \n',L22); +% % +% % end +%% +%{ +% printfig=1; +addpath Z:\Planar_Uncertainty_work\Results\final_plots\applyhatch_pluscolor_bundle\; +set(0,'DefaultAxesFontName', 'Times New Roman'); +crimson=rgb('crimson'); +dblue=rgb('deepskyblue'); +blueviolet=rgb('blueviolet'); +lw=1.5;fs=20;lw2=2.5; +% Plot uncertainty histograms +figure;hold on; +set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs); +bh=bar([covg1; covg2; covg3; covg4; covg5]); +set(bh(1),'Facecolor',crimson,'Edgecolor','k','LineWidth',lw2);set(bh(2),'Facecolor',crimson,'Edgecolor','k','LineStyle','--','LineWidth',lw2); +set(bh(3),'Facecolor',dblue,'Edgecolor','k','LineWidth',lw2);set(bh(4),'Facecolor',dblue,'Edgecolor','k','LineStyle','--','LineWidth',lw2); +set(bh(5),'Facecolor',blueviolet,'Edgecolor','k','LineWidth',lw2);set(bh(6),'Facecolor',blueviolet,'Edgecolor','k','LineStyle','--','LineWidth',lw2); +xl=xlim(gca); +plot(linspace(xl(1),xl(2),5),68.5.*ones(5,1),'k--'); +hold off; +set(gca,'XTick',[1 2 3 4 5]) +set(gca,'XTickLabel',casename_org); +box on; +set(gcf,'color','white'); +set(gcf,'units','normalized','outerposition',[0 0 0.8 0.8]); +ylabel('Coverage (%)'); +title('Standard Uncertainty Coverage'); +legend('MC_x','MC_y','PD_x','PD_y','CS_x','CS_y','orientation','horizontal'); +% applyhatch_pluscolor(gcf,'./././././././././././././././',1,360); +if printfig==1 + export_fig(gcf,fullfile(outputdir{2},'coverage_plot.png'),'-painters','-r360'); +end + +% Plot the uncertainty histograms +% plot_MC_uncertainty(A,casename,outputdir,printfig,runCS,runIM,runMC); +%} + +%Spatial uncertainties + + +%Uncertainty Timeseries + + + + +% end \ No newline at end of file diff --git a/postprocess_codes/socdiff.m b/postprocess_codes/socdiff.m new file mode 100644 index 0000000..f05c359 --- /dev/null +++ b/postprocess_codes/socdiff.m @@ -0,0 +1,56 @@ +function [N]=socdiff(N,dx,dir) + +size_N = size(N); + +if isempty(dir) + if size_N(1)~=1 + dir = 1; + elseif size_N(2)~=1 + dir = 2; + elseif size_N(3)~=1 + dir = 3; + else + dir = 1; + end +end + +ndim_N = ndims(N); +if ndim_N > 3 + error('function only defined up to 3D matrices') +end + +if isempty(dx) + dx = 1; +end + + +N = permute(N, [dir, 1:dir-1, dir+1:ndim_N]); +size_N = size(N); %find again for permuted matrix +ndim_N = ndims(N); + + +NI = size_N(1); +NJ = size_N(2); +if ndim_N == 3 + NK = size_N(3); +else + NK = 1; +end + + +if NI < 3 + error('stencil size is 3 pts in active direction') +end + +DN = zeros(NI,NJ,NK); +%for j=1:NJ +% for k=1:NK + DN(1 ,:,:) = ( -N( 3 ,:,:) +4*N( 2 ,:,:) -3*N(1 ,:,:)) / (2*dx); + DN(2:NI-1,:,:) = ( N( 3:NI,:,:) - N(1:NI-2,:,:)) / (2*dx); + DN(NI ,:,:) = (3*N(NI ,:,:) -4*N(NI-1,:,:) + N( NI-2,:,:)) / (2*dx); +% end +%end + +%keyboard + +N = ipermute(DN, [dir, 1:dir-1, dir+1:ndim_N]); diff --git a/postprocess_codes/stagnation_flow_postprocess.m b/postprocess_codes/stagnation_flow_postprocess.m new file mode 100644 index 0000000..23b1ba9 --- /dev/null +++ b/postprocess_codes/stagnation_flow_postprocess.m @@ -0,0 +1,399 @@ +% Testing IXX for STagnation Flow +clear all; +%Load true solution:- +Ut=load('Z:\Planar_Uncertainty_work\codes\matfiles\Utrue.mat'); +Vt=load('Z:\Planar_Uncertainty_work\codes\matfiles\Vtrue.mat'); +Utrue=Ut.Ue; +Vtrue=Vt.Ve; + + + +Un_fit=load('Z:\Planar_Uncertainty_work\codes\matfiles\fit_un.mat'); +unfitx=Un_fit.u_pUPe1; +unfity=Un_fit.u_pVPe1; +% figure;imagesc(Ut1);colorbar;caxis([-0.06 0.06]); +% figure;imagesc(Vt1);colorbar;caxis([0 3.5]); + +istring1=sprintf(['%%s%%0%0.0fd.','mat'],6); + +% dirsol='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\stagnation_flow\PIV_scc_finalstag_pass4_'; +% dirsol2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\stagnation_flow\PIV_scc_finalstag_pass5_'; + +% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\PIV_scc_finalstag_pass4_'; +% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\newixx\PIV_scc_finalstag_pass4_'; +% dirsol2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\PIV_scc_finalstag_pass5_'; +dirsol='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\stagnation_flow\PIV_scc_finalstag_pass4_'; +dirsol2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\stagnation_flow\ws32IM\PIV_scc_finalstag_pass4_'; + + +% outdir1='Z:\Planar_Uncertainty_work\Results\Ixx_with_cclsq_dia\plots\Stagnation_flow\'; +% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\plots\stagnation_flow\'; +% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\With_zero_center_gauss\plots\stagnation_flow\'; +% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\stagnation_flow\'; +% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\stagnation_flow\newixx\'; +outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\stagnation_flow\ws32IM\'; + +biasfact=0; + +if ~exist(outdir1,'dir') + mkdir(outdir1); +end +casename1='diacclsq2'; +outdir=[outdir1,casename1]; +printfig=1; + +t=99; + +Sx=64; +Sy=80; + +for i=1:t + sccsol=load(sprintf(istring1,dirsol,i)); + sccsol2=load(sprintf(istring1,dirsol2,i)); + + %Solution + Us(:,:,i)=sccsol.U; + Vs(:,:,i)=sccsol.V; + + %Error + err_sccu(:,:,i)=abs(Utrue-Us(:,:,i)); + err_sccv(:,:,i)=abs(Vtrue-Vs(:,:,i)); + + % Image Matching + Imx=reshape(sccsol2.uncertainty(:,9),Sx,Sy); + Imy=reshape(sccsol2.uncertainty(:,10),Sx,Sy); + + MI(:,:,i)=reshape(sccsol.uncertainty(:,5),Sx,Sy); + autod=reshape(sccsol.uncertainty(:,6),Sx,Sy); + Ixx=reshape(sccsol.uncertainty(:,7),Sx,Sy); + Iyy=reshape(sccsol.uncertainty(:,8),Sx,Sy); + + % correct gradient (if desired) + [Udiff]=socdiff(Us(:,:,i),16,1); + [Vdiff]=socdiff(Vs(:,:,i),16,2); + Ixx = real(sqrt(Ixx.^2 - (autod.^2/16).*(Udiff).^2 )); + Iyy = real(sqrt(Iyy.^2 - (autod.^2/16).*(Vdiff).^2 )); + + %Gcc peak location + Gpx=reshape(sccsol.uncertainty(:,15),Sx,Sy); + Gpy=reshape(sccsol.uncertainty(:,16),Sx,Sy); + + %Scc peak location + Spx=biasfact*reshape(sccsol.uncertainty(:,17),Sx,Sy); + Spy=biasfact*reshape(sccsol.uncertainty(:,18),Sx,Sy); + + %Bias error + mewx=abs(Spx-Gpx); + mewy=abs(Spy-Gpy); + + + %Scaling for ixx + Ixx=Ixx./sqrt(MI(:,:,i)); + Iyy=Iyy./sqrt(MI(:,:,i)); + + % Final second order moment uncertainty + IXX1(:,:,i)=sqrt(Ixx.^2+mewx.^2+unfitx.^2); + IYY1(:,:,i)=sqrt(Iyy.^2+mewy.^2+unfity.^2); + % Image Matching + IXX2(:,:,i)=sqrt(Imx.^2+unfitx.^2); + IYY2(:,:,i)=sqrt(Imy.^2+unfity.^2); + +end +%% convert to vectors +IXX1=IXX1(:); +IYY1=IYY1(:); +IXX2=IXX2(:); +IYY2=IYY2(:); +err_sccu=err_sccu(:); +err_sccv=err_sccv(:); +MI=MI(:); +%% +%zero value removal +inx=(find(IXX1(:)==0)); +iny=(find(IYY1(:)==0)); +%keyboard; +IXX1(inx)=[]; +IYY1(iny)=[]; +IXX2(inx)=[]; +IYY2(iny)=[]; +err_sccu(inx)=[]; +err_sccv(iny)=[]; +%% +%nan value removal +inx2=(find(isnan(IXX2))); +iny2=(find(isnan(IXX2))); + +%keyboard; +IXX1(inx2)=[]; +IYY1(iny2)=[]; +IXX2(inx2)=[]; +IYY2(iny2)=[]; +err_sccu(inx2)=[]; +err_sccv(iny2)=[]; + +inx3=(find(isnan(IXX1))); +iny3=(find(isnan(IXX1))); +IXX1(inx3)=[]; +IYY1(iny3)=[]; +IXX2(inx3)=[]; +IYY2(iny3)=[]; +err_sccu(inx3)=[]; +err_sccv(iny3)=[]; + +[length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)] +%% +%invalid vector removal +inx1=(find(err_sccu(:)>1)); +iny1=(find(err_sccv(:)>1)); +%keyboard; +IXX1(inx1)=[]; +IYY1(iny1)=[]; +IXX2(inx1)=[]; +IYY2(iny1)=[]; +err_sccu(inx1)=[]; +err_sccv(iny1)=[]; + +%% + +gmx=(find(err_sccu(:)<=1)); +gmy=(find(err_sccv(:)<=1)); +ngmx=length(gmx); +ngmy=length(gmy); +% N=r*s*t; +% [(N-ngmx)*100/N (N-ngmy)*100/N]; + +cnt=0; +for j=1:length(IXX1) +if err_sccu(j)<=IXX1(j) +cnt=cnt+1; +end +end +covx=100*cnt/ngmx; + +cnt=0; +for j=1:length(IYY1) +if err_sccv(j)<=IYY1(j) +cnt=cnt+1; +end +end +covy=100*cnt/ngmy; + +cnt=0; +for j=1:length(IXX2) +if err_sccu(j)<=IXX2(j) +cnt=cnt+1; +end +end +covx2=100*cnt/ngmx; + +cnt=0; +for j=1:length(IYY2) +if err_sccv(j)<=IYY2(j) +cnt=cnt+1; +end +end +covy2=100*cnt/ngmy; + +[covx covx2]% +[covy covy2]% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%% +lw=3; +fs=18; + %% Plot error and uncertainty histogram +thresh=0.3; + +vec=linspace(0,thresh,60); + +Nrex=histc(err_sccu,vec); +Nrey=histc(err_sccv,vec); +Nrdx=histc(IXX1,vec); +Nrdy=histc(IYY1,vec); +Nrdx2=histc(IXX2,vec); +Nrdy2=histc(IYY2,vec); + + +ll=1:100:max(Nrex); +ll2=1:100:max(Nrey); + + +sigu=rms(err_sccu(err_sccu<=thresh)).*ones(length(ll),1); +sigv=rms(err_sccv(err_sccv<=thresh)).*ones(length(ll2),1); + +mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1); +miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + +mixx2=rms(IXX2(IXX2<=thresh)).*ones(length(ll),1); +miyy2=rms(IYY2(IYY2<=thresh)).*ones(length(ll2),1); + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--'); +hold off; + +% title('Error and Uncertainty histogram(X)'); +title(sprintf(['Ixx=',num2str(covx),' IMx=',num2str(covx2)])); +%title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 thresh 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outdir,'Uhist2.png']); +end + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--'); +hold off; + +% title('Error and Uncertainty histogram(Y)'); +title(sprintf(['Iyy=',num2str(covy),' IMy=',num2str(covy2)])); +%title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 thresh 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outdir,'Vhist2.png']); +end + +fprintf(['RMSerr X','\n']); +fprintf(num2str([sigu(1) mixx1(1) mixx2(1)],'%02.4f\t')); +fprintf(['\nRMSerr Y','\n']); +fprintf(num2str([sigv(1) miyy1(1) miyy2(1)],'%02.4f\t')); +fprintf('\n'); +%% +plotflag=1; + +if plotflag==1; + + Nbin=10; + + maxixx=0.15; + minixx=0.02; + vu=linspace(minixx,maxixx,Nbin); + maxiyy=0.3; + miniyy=0.04; + vv=linspace(miniyy,maxiyy,Nbin); + + rmserru1=zeros(1,length(vu)); + rmserrv1=zeros(1,length(vv)); + getixx1=zeros(1,length(vu)); + getiyy1=zeros(1,length(vv)); + numcx1=zeros(1,length(vu)); + numcy1=zeros(1,length(vv)); + covxbin1=zeros(1,length(vu)); + covybin1=zeros(1,length(vv)); + + rmserru2=zeros(1,length(vu)); + rmserrv2=zeros(1,length(vv)); + getixx2=zeros(1,length(vu)); + getiyy2=zeros(1,length(vv)); + numcx2=zeros(1,length(vu)); + numcy2=zeros(1,length(vv)); + covxbin2=zeros(1,length(vu)); + covybin2=zeros(1,length(vv)); + + %p=1; + + for p=1:length(vu) + if p=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1ax & xvecay & yvecthresh1)); +iny1=(find(error_V(:)>thresh1)); +%keyboard; +IXX1(inx1)=[]; +IYY1(iny1)=[]; +Imx(inx1)=[]; +Imy(iny1)=[]; +error_U(inx1)=[]; +error_V(iny1)=[]; +noabserrx(inx1)=[]; +noabserry(iny1)=[]; +%% +figure;hist(noabserrx,60); +figure;hist(noabserry,60); + +[mean(noabserrx) mean(noabserry)] +%% Find Coverage + +gmx=(find(error_U(:)<=1)); +gmy=(find(error_V(:)<=1)); +ngmx=length(gmx); +ngmy=length(gmy); + +%keyboard; +% N=r*s*t; +% [(N-ngmx)*100/N (N-ngmy)*100/N]; + +cnt=0; +for j=1:length(IXX1) +if error_U(j)<=IXX1(j) +cnt=cnt+1; +end +end +covx=100*cnt/ngmx; + +cnt=0; +for j=1:length(IYY1) +if error_V(j)<=IYY1(j) +cnt=cnt+1; +end +end +covy=100*cnt/ngmy; + +cnt=0; +for j=1:length(Imx) +if error_U(j)<=Imx(j) +cnt=cnt+1; +end +end +covx2=100*cnt/ngmx; + +cnt=0; +for j=1:length(Imy) +if error_V(j)<=Imy(j) +cnt=cnt+1; +end +end +covy2=100*cnt/ngmy; + +[covx covx2] +[covy covy2] + +% [rms(error_U) rms(error_V)] + +% IXX2=Imx; +% IYY2=Imy; + +% save('VRSOM.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2'); + + %% Plot error and uncertainty histogram +thresh=0.2; + +vec=linspace(0,thresh,60); + +Nrex=histc(error_U,vec); +Nrey=histc(error_V,vec); +Nrdx=histc(IXX1,vec); +Nrdy=histc(IYY1,vec); +Nrdx2=histc(Imx,vec); +Nrdy2=histc(Imy,vec); + + +ll=1:100:max(Nrex); +ll2=1:100:max(Nrey); + + +sigu=rms(error_U(error_U<=thresh)).*ones(length(ll),1); +sigv=rms(error_V(error_V<=thresh)).*ones(length(ll2),1); + +mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1); +miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1); + +mixx2=rms(Imx(Imx<=thresh)).*ones(length(ll),1); +miyy2=rms(Imy(Imy<=thresh)).*ones(length(ll2),1); + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--'); +hold off; + +% title('Error and Uncertainty histogram(X)'); +title(sprintf(['Ixx=',num2str(covx),' IMx=',num2str(covx2)])); +%title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 0.2 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outdir,'Uhist2.png']); +end + +figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on; +plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--'); +hold off; + +% title('Error and Uncertainty histogram(Y)'); +title(sprintf(['Iyy=',num2str(covy),' IMy=',num2str(covy2)])); +%title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)])); +ylabel('No. of Data Points');xlabel('pixels'); +axis([0 0.2 0 inf]); +if printfig==1 + print(gcf,'-dpng',[outdir,'Vhist2.png']); +end + +fprintf(['RMSerr X','\n']); +fprintf(num2str([sigu(1) mixx1(1) mixx2(1)],'%02.4f\t')); +fprintf(['\nRMSerr Y','\n']); +fprintf(num2str([sigv(1) miyy1(1) miyy2(1)],'%02.4f\t')); +fprintf('\n'); +%% +plotflag=1; + +if plotflag==1; + + IXX2=Imx; + IYY2=Imy; + err_sccu=error_U; + err_sccv=error_V; + + Nbin=15; + + maxixx=0.1; + minixx=0; + vu=linspace(minixx,maxixx,Nbin); + maxiyy=0.1; + miniyy=0; + vv=linspace(miniyy,maxiyy,Nbin); + + rmserru1=zeros(1,length(vu)); + rmserrv1=zeros(1,length(vv)); + getixx1=zeros(1,length(vu)); + getiyy1=zeros(1,length(vv)); + numcx1=zeros(1,length(vu)); + numcy1=zeros(1,length(vv)); + covxbin1=zeros(1,length(vu)); + covybin1=zeros(1,length(vv)); + + rmserru2=zeros(1,length(vu)); + rmserrv2=zeros(1,length(vv)); + getixx2=zeros(1,length(vu)); + getiyy2=zeros(1,length(vv)); + numcx2=zeros(1,length(vu)); + numcy2=zeros(1,length(vv)); + covxbin2=zeros(1,length(vu)); + covybin2=zeros(1,length(vv)); + + %p=1; + + for p=1:length(vu) + if p=vu(p))); + [val2]= find((IXX2(:)=vu(p))); + elseif p==length(vu) + [val1]= find((IXX1(:)>=vu(p))); + [val2]= find((IXX2(:)>=vu(p))); + end + numcx1(p)=(length(val1)); + numcx2(p)=(length(val2)); + tmp1=(err_sccu(val1)); + tmp2=(err_sccu(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempixx1=IXX1(val1); + tempixx2=IXX2(val2); + cvrx1=find(tmp1=vv(p))); + [val2]= find((IYY2(:)=vv(p))); + elseif p==length(vv) + [val1]= find((IYY1(:)>=vv(p))); + [val2]= find((IYY2(:)>=vv(p))); + end + numcy1(p)=(length(val1)); + numcy2(p)=(length(val2)); + tmp1=(err_sccv(val1)); + tmp2=(err_sccv(val2)); + % tmp1(tmp1>1)=[]; + % tmp2(tmp2>1)=[]; + tempiyy1=IYY1(val1); + tempiyy2=IYY2(val2); + cvry1=find(tmp1