diff --git a/vortex_ring_postprocess_codes/Postprocess_VR_davis_01_11_2016.m b/vortex_ring_postprocess_codes/Postprocess_VR_davis_01_11_2016.m
new file mode 100644
index 0000000..9de7cd1
--- /dev/null
+++ b/vortex_ring_postprocess_codes/Postprocess_VR_davis_01_11_2016.m
@@ -0,0 +1,731 @@
+%Postprocess_stereo_uncertainty_vortexring thicksheet
+
+%Uncertainty Propagation for Stereo Geometric Reconstruction(Willert
+%Method)for vortex ring
+
+clear all;
+lw=4;fs=18;
+if ispc
+    main_dir='Z:\Stereo_Uncertainty_Work\';
+    addpath Z:\Stereo_Uncertainty_Work\Codes\New_codes\;
+
+else
+    main_dir='/home/shannon/a/bhattac3/Stereo_Uncertainty_Work/';
+    addpath /home/shannon/a/bhattac3/Stereo_Uncertainty_Work/Codes/New_codes/;
+
+end
+
+meth='CS';
+
+printfig=0;
+ensdispfldflag=0;
+
+wres=32;
+
+cam1='1';cam2='3';
+% cam1='2';cam2='4';
+
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    casename='CScam13';
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    casename='CScam24';
+end
+
+
+dirout=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','plots_final','Wihselfcal_angle_Davis');
+if ~exist(dirout,'dir')
+    mkdir(dirout);
+end
+dirout=[dirout,filesep,casename,'_'];
+
+
+stdjob= load(fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','jobs',['cam',cam1,cam2,'_VR_prana_fulljob_withselfcal.mat']));
+caljob=stdjob.caljobfile.datasave.caljob;
+calmat=[caljob.aXcam1 caljob.aYcam1 caljob.aXcam2 caljob.aYcam2];
+calmat;
+planarjob=stdjob.pranajob.Data;
+
+dispfield=load(fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Willert','dispfield_ensemble',['VR',cam1,cam2,'_scc_pass2_00025.mat']));
+
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    imagedir.imdir1=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images1','E_camera_1_frame_');
+    imagedir.imdir2=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images2','E_camera_3_frame_');
+
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    imagedir.imdir1=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images1','Camera_02_frame_');
+    imagedir.imdir2=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images2','Camera_04_frame_');
+
+end
+
+% Load Davis fields
+Davis13= load(fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Davis_reconstruction','Vortex_Ring',['cam',cam1,cam2],'Willert',['cam',cam1,'cam',cam2,'_willert_recon_field.mat']));
+% Davis13= load('Z:\Stereo_Uncertainty_Work\New_Stereo_uncertainty_tests\Vortex_Ring\Davis_reconstruction\Vortex_Ring\cam24\Willert\cam2cam4_willert_recon_field.mat');
+
+
+%% Input Directories
+basedir=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Willert',filesep);
+
+dirwc1=[basedir,'dispfields',filesep,'VR',cam1,cam2,'_scc_pass4_'];
+istring1=sprintf(['%%s%%0%0.0fd.','mat'],5);
+istring2=sprintf(['%%s%%0%0.0fd.','mat'],5);
+
+%Camera1 directory
+dircam1a=[basedir,'Camera',cam1,filesep,'VR',cam1,cam2,'_scc_pass4_'];
+
+%Camera2 directory
+dircam2a=[basedir,'Camera',cam2,filesep,'VR',cam1,cam2,'_scc_pass4_'];
+
+%reconstructed field directories
+dirrec=[basedir,'Camera',cam1,'Camera',cam2,'_3Cfields',filesep,'piv_2d3c_cam',cam1,'cam',cam2,'_pass_4_'];
+
+%true solution directory
+dirtrue=fullfile(main_dir,'VortexRing_True_solution','MATLAB_Data','Eval_1','E_eval1_TrueSolution_');
+%dirtrue='/home/shannon/a/bhattac3/Stereo_Uncertainty_Work/Process_vortexring_cam2cam4/Velocity_field_with_selfcal/Willert/Camera2Camera4_3Cfields/piv_2d3c_cam2cam4_pass_4_';
+
+scaling=stdjob.scaling.wil;
+mx=scaling.xscale;
+my=scaling.yscale;
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    mz=scaling.yscale;
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    mz=scaling.xscale;
+end
+%seconds per frame
+spf=0.001;
+
+t=50;   
+%get uncertainties in x and y
+ofs=24;
+%% Find uncertainty in angles  
+outputdir=dirout;
+fstart=ofs;
+Nim=t;
+
+[Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,alpha1,alpha2,beta1,beta2,Umx,Umy]=stereo_calcoeff_uncertainty(caljob,planarjob,dispfield,imagedir,outputdir,ensdispfldflag,wres,fstart,Nim);
+ 
+% [mean(alpha1(:)) mean(alpha2(:)) mean(beta1(:)) mean(beta2(:))]
+% 
+% [mean(Un_alpha1(:)) mean(Un_alpha2(:)) mean(Un_beta1(:)) mean(Un_beta2(:))]
+
+%% 
+[Sx,Sy]=size(alpha1);
+
+for i=1:t
+  
+   %% Load reconstructed solution
+   Us(:,:,i)=(spf/mx)*1e3*Davis13.Uw(1:end-1,1:end-1,i);
+   Vs(:,:,i)=(spf/my)*1e3*Davis13.Vw(1:end-1,1:end-1,i);
+   Ws(:,:,i)=(spf/mz)*1e3*Davis13.Ww(1:end-1,1:end-1,i);
+
+   xgrid=1000*Davis13.Xw(1:end-1,1:end-1);
+   ygrid=1000*Davis13.Yw(1:end-1,1:end-1);
+   zgrid=1000*Davis13.Zw(1:end-1,1:end-1);
+   
+   [r,c]=size(xgrid);
+   
+   %% Load true solution
+
+   TS1=load(sprintf(istring2,dirtrue,i+24));
+   TS=TS1.velocity_data;
+   ut=(spf/mx)*1000*TS.x_velocity;
+   vt=(spf/my)*1000*TS.y_velocity;
+   wt=(spf/mz)*1000*TS.z_velocity;
+   [xt,yt]=meshgrid(-34:0.45:20,-27:0.45:27);
+%    xt=TS.x_coordinate';
+%    yt=TS.y_coordinate';
+%    zt=TS.z_coordinate';
+%keyboard;
+   
+   Ut(:,:,i)=interp2(xt,yt,ut,xgrid,ygrid,'spline',0);
+   Vt(:,:,i)=interp2(xt,yt,vt,xgrid,ygrid,'spline',0);
+   Wt(:,:,i)=interp2(xt,yt,wt,xgrid,ygrid,'spline',0);
+
+%keyboard;
+   
+%    figure(4);
+%    subplot(2,3,1);imagesc(xgrid(:),ygrid(:),Us(:,:,i));colorbar;title('Us');axis equal square;caxis([-1 3]);
+%    subplot(2,3,2);imagesc(xgrid(:),ygrid(:),Vs(:,:,i));colorbar;title('Vs');axis equal square;caxis([-1.5 1.5]);
+%    subplot(2,3,3);imagesc(xgrid(:),ygrid(:),Ws(:,:,i));colorbar;title('Ws');axis equal square;caxis([-1 1]);
+%    subplot(2,3,4);imagesc(xgrid(:),ygrid(:),Ut(:,:,i));colorbar;title('Ut');axis equal square;caxis([-1 3]);
+%    subplot(2,3,5);imagesc(xgrid(:),ygrid(:),Vt(:,:,i));colorbar;title('Vt');axis equal square;caxis([-1.5 1.5]);
+%    subplot(2,3,6);imagesc(xgrid(:),ygrid(:),Wt(:,:,i));colorbar;title('Wt');axis equal square;caxis([-1 1]);
+%    keyboard;
+% %    pause(0.05);
+   
+   %% Get error fields
+   errx1(:,:,i)=(Ut(:,:,i)-Us(:,:,i));
+   erry1(:,:,i)=(Vt(:,:,i)-Vs(:,:,i));
+   errz1(:,:,i)=(Wt(:,:,i)-Ws(:,:,i));
+   
+   uttemp=Ut(:,:,i);
+   vttemp=Vt(:,:,i);
+   ustemp=Us(:,:,i);
+   vstemp=Vs(:,:,i);
+   
+   truediv(:,:,i)=divergence(xgrid,ygrid,uttemp,vttemp);
+   soldiv(:,:,i)=divergence(xgrid,ygrid,ustemp,vstemp);
+   %% Load cam1cam2 2d fields
+    %get 2d velocities
+   U1=Davis13.Us1(1:end-1,1:end-1,i);
+   V1=Davis13.Vs1(1:end-1,1:end-1,i);
+   X1=Davis13.X(1:end-1,1:end-1);
+   Y1=Davis13.Y(1:end-1,1:end-1);
+   
+   U2=Davis13.Us2(1:end-1,1:end-1,i);
+   V2=Davis13.Vs2(1:end-1,1:end-1,i);
+   X2=Davis13.X(1:end-1,1:end-1);
+   Y2=Davis13.Y(1:end-1,1:end-1);
+   
+   %Get 2d uncertainties
+   %Camera1
+   imx1(:,:,i)=Davis13.UCSx1(1:end-1,1:end-1,i);%;
+   imy1(:,:,i)=Davis13.UCSy1(1:end-1,1:end-1,i);%
+   
+   %Camera2
+   imx2(:,:,i)=Davis13.UCSx2(1:end-1,1:end-1,i);%
+   imy2(:,:,i)=Davis13.UCSy2(1:end-1,1:end-1,i);%
+  
+  
+   %% Stereo Uncertainty Propagation
+
+    Unu1=imx1(:,:,i);
+    Unu2=imx2(:,:,i);
+    Unv1=imy1(:,:,i);
+    Unv2=imy2(:,:,i);
+
+   [Un_u1,Un_v1,Un_w1,JU,JV,JW]=stereo_uncertainty_propagation_withcalcoeff(Unu1,Unv1,Unu2,Unv2,U1,V1,U2,V2,Ws(:,:,i),Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,alpha1,alpha2,beta1,beta2,mx,my);
+
+   [Un_u3,Un_v3,Un_w3]=stereo_uncertainty_propagation_noangle(Unu1,Unv1,Unu2,Unv2,alpha1,alpha2,beta1,beta2,mx,my);
+   
+   % Uncertainty in magnification
+%    Umagx(:,:,i)=sqrt((Umx^2)*((Us(:,:,i)*mx).^2)/(mx^4));
+%    Umagy(:,:,i)=sqrt((Umy^2)*((Vs(:,:,i)*my).^2)/(my^4));
+%    Umagz(:,:,i)=sqrt((Umx^2)*((Ws(:,:,i)*mx).^2)/(mx^4));
+   
+%     unim_u(:,:,i)=sqrt(Un_u1.^2+Umagx(:,:,i).^2);
+%     unim_v(:,:,i)=sqrt(Un_v1.^2+Umagy(:,:,i).^2);
+%     unim_w(:,:,i)=sqrt(Un_w1.^2+Umagz(:,:,i).^2);
+%     unim_um(:,:,i)=Un_u3;
+%     unim_vm(:,:,i)=Un_v3;
+%     unim_wm(:,:,i)=Un_w3;
+   
+    unim_u(:,:,i)=Un_u1;
+    unim_v(:,:,i)=Un_v1;
+    unim_w(:,:,i)=Un_w1;
+    unim_um(:,:,i)=Un_u3;
+    unim_vm(:,:,i)=Un_v3;
+    unim_wm(:,:,i)=Un_w3;
+
+%     unim_u(:,:,i)=Un_u3;
+%     unim_v(:,:,i)=Un_v3;
+%     unim_w(:,:,i)=Un_w3;
+%     unim_um(:,:,i)=Un_u1;
+%     unim_vm(:,:,i)=Un_v1;
+%     unim_wm(:,:,i)=Un_w1;
+
+    %% RMS error and uncertainty time series
+%     tex1=errx1(:,:,i);tux1=Un_u1(:);
+%     tex1=tex1(:);    nx=find(abs(tex1)>0.5);
+%     tex1(nx)=[];
+%     tux1(nx)=[];
+%     siguerr(i)=rms(tex1);siguun(i)=rms(tux1);
+%     
+%     tey1=erry1(:,:,i);tuy1=Un_v1(:);
+%     tey1=tey1(:);    ny=find(abs(tey1)>0.5);
+%     tey1(ny)=[];
+%     tuy1(ny)=[];
+%     sigverr(i)=rms(tey1);sigvun(i)=rms(tuy1);
+%     
+%     tez1=errz1(:,:,i);tuz1=Un_w1(:);
+%     tez1=tez1(:);    nz=find(abs(tez1)>2);
+%     tez1(nz)=[];
+%     tuz1(nz)=[];
+%     sigwerr(i)=rms(tez1);sigwun(i)=rms(tuz1);
+
+end
+soldiv2=soldiv(abs(soldiv)<1);
+[mean(truediv(:)) mean(soldiv2(:))]
+[rms(truediv(:)) rms(soldiv2(:))]
+keyboard;
+%Good measurement cut off
+ulx=0.5;uly=0.5;ulz=1.5;
+%% Plot RMS Uncertainty and Error fields (Spatial RMS Profile)
+
+for i=1:Sx
+    for j=1:Sy
+        
+        
+        rms2derrx=squeeze(abs(errx1(i,j,:)));
+        rms2dunx=squeeze(abs(unim_u(i,j,:)));
+        ix1=find(rms2derrx>ulx);
+        rms2derrx(ix1)=[];
+        rms2dunx(ix1)=[];
+%         Nx(i,j)=length(ix1);
+        
+        rms2derry=squeeze(abs(erry1(i,j,:)));
+        rms2duny=squeeze(abs(unim_v(i,j,:)));
+        iy1=find(rms2derry>uly);
+        rms2derry(iy1)=[];
+        rms2duny(iy1)=[];
+%         Ny(i,j)=length(iy1);
+        
+        rms2derrz=squeeze(abs(errz1(i,j,:)));
+        rms2dunz=squeeze(abs(unim_w(i,j,:)));
+        iz1=find(rms2derrz>ulz);
+        rms2derrz(iz1)=[];
+        rms2dunz(iz1)=[];
+%         Nz(i,j)=length(iz1);
+       
+        
+        Ue(i,j)=rms(rms2derrx);
+        Uun(i,j)=rms(rms2dunx);
+        Ve(i,j)=rms(rms2derry);
+        Vun(i,j)=rms(rms2duny);
+        We(i,j)=rms(rms2derrz);
+        Wun(i,j)=rms(rms2dunz);
+        
+    end
+end
+
+figure;
+subplot(2,3,1);imagesc(Ue);caxis([0 0.2]);
+subplot(2,3,2);imagesc(Ve);caxis([0 0.2]);
+subplot(2,3,3);imagesc(We);caxis([0 0.3]);
+subplot(2,3,4);imagesc(Uun);caxis([0 0.2]);
+subplot(2,3,5);imagesc(Vun);caxis([0 0.2]);
+subplot(2,3,6);imagesc(Wun);caxis([0 0.3]);
+%% Plot RMS Uncertainty and Error fields (Temporal RMS Profile)
+
+for k=1:t
+    
+        rex=squeeze(abs(errx1(:,:,k)));
+        rIMx=squeeze(abs(unim_u(:,:,k)));
+        rex=rex(:);
+        rIMx=rIMx(:);
+        ix1=find(rex>ulx);
+        rex(ix1)=[];
+        rIMx(ix1)=[];
+%         Nx(i,j)=length(ix1);
+        
+        rey=squeeze(abs(erry1(:,:,k)));
+        rIMy=squeeze(abs(unim_v(:,:,k)));
+        rey=rey(:);
+        rIMy=rIMy(:);
+        iy1=find(rey>uly);
+        rey(iy1)=[];
+        rIMy(iy1)=[];
+%         Ny(i,j)=length(iy1);
+        
+        rez=squeeze(abs(errz1(:,:,k)));
+        rIMz=squeeze(abs(unim_w(:,:,k)));
+        rez=rez(:);
+        rIMz=rIMz(:);
+        iz1=find(rez>ulz);
+        rez(iz1)=[];
+        rIMz(iz1)=[];
+%         Nz(i,j)=length(iz1);
+       
+        
+        siguerr(k)=rms(rex);
+        siguun(k)=rms(rIMx);
+        sigverr(k)=rms(rey);
+        sigvun(k)=rms(rIMy);
+        sigwerr(k)=rms(rez);
+        sigwun(k)=rms(rIMz);
+        
+   
+end
+
+figure;hold on;
+plot(1:t,siguerr,'r-o',1:t,siguun,'r-*');
+plot(1:t,sigverr,'b-o',1:t,sigvun,'b-*');
+plot(1:t,sigwerr,'g-o',1:t,sigwun,'g-*');
+hold off;
+title('RMS Time Series (Cam24)');xlabel('Frame No.');
+ylabel('RMS error, uncertainty(pix/frame)');
+% hleg=legend({'\sigma^{e}_{u}','\sigma^{CS}_{u}','\sigma^{e}_{v}','\sigma^{CS}_{v}','\sigma^{e}_{w}','\sigma^{CS}_{w}'});
+% set(hleg,'location','eastoutside');
+set(gca,'FontSize',20);
+axis([0 50 0 0.5]);axis square;
+%% RMS 2D Uncertainty spatial maps
+nc=20;
+IMX1=rms(imx1,3);
+IMY1=rms(imy1,3);
+IMX2=rms(imx2,3);
+IMY2=rms(imy2,3);
+figure;
+subplot(2,2,1);contour(xgrid,ygrid,IMX1,nc);caxis([0 0.3]);title('\sigma^{IM}_{U1}');axis square;
+subplot(2,2,2);contour(xgrid,ygrid,IMY1,nc);caxis([0 0.3]);title('\sigma^{IM}_{V1}');axis square;
+subplot(2,2,3);contour(xgrid,ygrid,IMX2,nc);caxis([0 0.3]);title('\sigma^{IM}_{U2}');axis square;
+subplot(2,2,4);contour(xgrid,ygrid,IMY2,nc);caxis([0 0.3]);title('\sigma^{IM}_{V2}');axis square;
+
+%%
+save([meth,'_Cam',cam1,cam2,'_RMSprofiles.mat'],'xgrid','ygrid','IMX1',...
+    'IMX2','IMY1','IMY2','siguerr','siguun','sigverr','sigvun','sigwerr','sigwun',...
+    'Ue','Uun','Ve','Vun','We','Wun');
+
+keyboard;
+
+%% plot rms time series
+% figure;hold on;
+% plot(1:t,siguerr,'r-o',1:t,siguun,'r-*');
+% plot(1:t,sigverr,'b-o',1:t,sigvun,'b-*');
+% plot(1:t,sigwerr,'g-o',1:t,sigwun,'g-*');
+% hold off;
+% title('RMS Time Series (Cam24)');xlabel('Frame No.');
+% ylabel('RMS error, uncertainty(pix/frame)');
+% % hleg=legend({'\sigma^{e}_{u}','\sigma^{CS}_{u}','\sigma^{e}_{v}','\sigma^{CS}_{v}','\sigma^{e}_{w}','\sigma^{CS}_{w}'});
+% % set(hleg,'location','eastoutside');
+% set(gca,'FontSize',20);
+% axis([0 50 0 0.5]);axis square;
+% print(gcf,'-dpng','Cam24CSeu','-r300');
+% 
+% keyboard;
+
+%% Plot RMS Uncertainty and Error fields
+%{
+Uun=rms(unim_u,3);
+Vun=rms(unim_v,3);
+Wun=rms(unim_w,3);
+
+Ue=rms(abs(errx1),3);
+Ve=rms(abs(erry1),3);
+We=rms(abs(errz1),3);
+
+figure;
+subplot(2,3,1);imagesc(Ue);caxis([0 0.2]);
+subplot(2,3,2);imagesc(Ve);caxis([0 0.2]);
+subplot(2,3,3);imagesc(We);caxis([0 0.3]);
+subplot(2,3,4);imagesc(Uun);caxis([0 0.2]);
+subplot(2,3,5);imagesc(Vun);caxis([0 0.2]);
+subplot(2,3,6);imagesc(Wun);caxis([0 0.3]);
+%}
+% keyboard;
+%%
+%[mean(errx(:)) mean(erry(:)) mean(errz(:))]
+
+% [mean(imx1(:)) mean(imy1(:)) mean(Nim1(:)) mean(imx2(:)) mean(imy2(:)) mean(Nim2(:))]
+% [rms(imx1(:)) rms(imy1(:)) mean(Nim1(:)) rms(imx2(:)) rms(imy2(:)) mean(Nim2(:))]
+
+%%
+errx=(errx1);
+erry=(erry1);
+errz=(errz1);
+
+% errx1(abs(errx1)>0.3)=[];
+% erry1(abs(erry1)>0.3)=[];
+% errz1(abs(errz1)>0.3)=[];
+% 
+% Nbin1=60;
+% figure;hist(errx1(:),Nbin1);
+% figure;hist(erry1(:),Nbin1);
+% figure;hist(errz1(:),Nbin1);
+
+Nc=100;
+
+ul=0.5;
+
+vx1=linspace(0,ul,Nc);
+vy1=linspace(0,ul,Nc);
+
+vx2=linspace(0,ul,Nc);
+vy2=linspace(0,ul,Nc);
+
+Nu1=histc(imx1(:),vx1);
+Nv1=histc(imy1(:),vy1);
+Nu2=histc(imx2(:),vx2);
+Nv2=histc(imy2(:),vy2);
+
+figure;set(gcf,'DefaultLineLineWidth',1);set(gca,'FontSize',fs);
+plot(vx1,Nu1,'r-o',vy1,Nv1,'b-o',vx2,Nu2,'r-+',vy2,Nv2,'b-+');
+title('Planar Correlation Uncertainties');
+legend('U1','V1','U2','V2');
+if printfig==1
+    print(gcf,'-dpng',[dirout,'planarhist.png'],'-r300');
+end
+%% RMS error and uncertainty field
+
+% % rex=rms(errx,3);
+% % rey=rms(erry,3);
+% % rez=rms(errz,3);
+% % rIMx=rms(unim_u,3);
+% % rIMy=rms(unim_v,3);
+% % rIMz=rms(unim_w,3);
+% % 
+% % figure;
+% % imagesc(rex(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{e}_{u}');axis square tight;
+% % set(gca,'FontSize',20);
+% % print(gcf,'-dpng','sigerru.png','-r300');
+% % figure;
+% % imagesc(rIMx(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{CS}_{u}');axis square tight;
+% % set(gca,'FontSize',20);
+% % print(gcf,'-dpng','sigUnu.png','-r300');
+% % 
+% % figure;imagesc(rez(end:-1:1,:));caxis([0 0.5]);axis off;title('\sigma^{e}_{w}');axis square tight;
+% % figure;imagesc(rIMz(end:-1:1,:));caxis([0 0.5]);axis off;title('\sigma^{CS}_{w}');axis square tight;
+% % figure;imagesc(rey(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{e}_{v}');axis square tight;
+% % figure;imagesc(rIMy(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{CS}_{v}');axis square tight;
+% % keyboard;
+%% eliminate invalid measurements
+errx=errx(:);
+erry=erry(:);
+errz=errz(:);
+unim_u=unim_u(:);
+unim_v=unim_v(:);
+unim_w=unim_w(:);
+
+
+unim_um=unim_um(:);
+unim_vm=unim_vm(:);
+unim_wm=unim_wm(:);
+
+% %Good measurement cut off
+% ulx=0.5;uly=0.5;ulz=1.5;
+
+inx=find(abs(errx(:))>ulx);
+iny=find(abs(erry(:))>uly);
+inz=find(abs(errz(:))>ulz);
+
+errx(inx)=[];
+erry(iny)=[];
+errz(inz)=[];
+unim_u(inx)=[];
+unim_v(iny)=[];
+unim_w(inz)=[];
+
+
+unim_um(inx)=[];
+unim_vm(iny)=[];
+unim_wm(inz)=[];
+
+% [mean(errx(:)) mean(erry(:)) mean(errz(:))]
+%%
+
+% Calculate coverage
+cnt1=0;cnt2=0;cnt3=0;cnt4=0;cnt5=0;cnt6=0;
+l1=length(errx);
+l2=length(erry);
+l3=length(errz);
+for m1=1:l1
+    if abs(errx(m1))<=unim_u(m1)
+        cnt1=cnt1+1;
+    end
+    if abs(errx(m1))<=unim_um(m1)
+        cnt4=cnt4+1;
+    end
+    
+end
+
+for m2=1:l2
+    if abs(erry(m2))<=unim_v(m2)
+        cnt2=cnt2+1;
+    end
+    if abs(erry(m2))<=unim_vm(m2)
+        cnt5=cnt5+1;
+    end
+end
+
+for m3=1:l3
+    
+    if abs(errz(m3))<=unim_w(m3)
+        cnt3=cnt3+1;
+    end
+    if abs(errz(m3))<=unim_wm(m3)
+        cnt6=cnt6+1;
+    end
+    
+end
+covx1=100*cnt1/l1;
+covy1=100*cnt2/l2;
+covz1=100*cnt3/l3;
+fprintf('\n Coverage with angle uncertainty');
+[covx1 covy1 covz1]
+covx2=100*cnt4/l1;
+covy2=100*cnt5/l2;
+covz2=100*cnt6/l3;
+fprintf('\n Coverage without angle uncertainty');
+[covx2 covy2 covz2]
+
+%%
+%% Plot Histograms
+if ispc
+addpath Z:\Planar_Uncertainty_work\codes\export_fig\;
+addpath Z:\Planar_Uncertainty_work\codes\color_code\;
+else
+    addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/export_fig/;
+    addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/color_code/;
+end
+
+set(0,'DefaultAxesFontName', 'Times New Roman');
+
+%Plot histograms
+Nb=30;
+
+vecx=linspace(0,ulx,Nb);
+vecy=linspace(0,uly,Nb);
+vecz=linspace(0,ulz,Nb);
+vecx1=linspace(-ulx,ulx,2*Nb-1);
+vecy1=linspace(-uly,uly,2*Nb-1);
+vecz1=linspace(-ulz,ulz,2*Nb-1);
+
+Nex=histc(errx,vecx1)/length(errx);
+Nimx=histc(unim_u,vecx)/length(errx);
+
+Ney=histc(erry,vecy1)/length(erry);
+Nimy=histc(unim_v,vecy)/length(erry);
+
+Nez=histc(errz,vecz1)/length(errz);
+Nimz=histc(unim_w,vecz)/length(errz);
+
+
+% X uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecx1,Nex,'-k');
+plot(vecx,Nimx,'r-');
+% HL1=fill(vecx1,Nex,'-k');
+% HL2=fill(vecx,Nimx,'r-');
+% set(HL1,'facealpha',.5);
+% set(HL2,'facealpha',.5);
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigu=rms(errx(abs(errx)<=ulx)).*ones(length(ll),1);
+mixx1=rms(unim_u(unim_u<=ulx)).*ones(length(ll),1);
+plot(sigu,ll,'k--',mixx1,ll,'r--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{x}',['\sigma_{x}^{',meth,'}']);
+axis([-ulx ulx 0 inf]);
+% title(['U Uncertainty Histogram coverage= ',num2str(covx1,'%02.02f')]);
+title('U Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Uhist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Uhist.png']),'-painters','-r360');
+end
+
+% Y uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecy1,Ney,'-k');
+plot(vecy,Nimy,'b-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigv=rms(erry(abs(erry)<=uly)).*ones(length(ll),1);
+miyy1=rms(unim_v(unim_v<=uly)).*ones(length(ll),1);
+plot(sigv,ll,'k--',miyy1,ll,'b--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{y}',['\sigma_{y}^{',meth,'}']);
+axis([-uly uly 0 inf]);
+% title(['V uncertainty histogram coverage= ',num2str(covy1,'%02.02f')]);
+title('V Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Vhist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Vhist.png']),'-painters','-r360');
+end
+
+% Z uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecz1,Nez,'-k');
+plot(vecz,Nimz,'g-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigw=rms(errz(abs(errz)<=ulz)).*ones(length(ll),1);
+mizz1=rms(unim_w(unim_w<=ulz)).*ones(length(ll),1);
+plot(sigw,ll,'k--',mizz1,ll,'g--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{z}',['\sigma_{z}^{',meth,'}']);
+axis([-ulz ulz 0 inf]);
+% title(['W uncertainty histogram coverage= ',num2str(covz1,'%02.02f')]);
+title('W Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Whist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Whist.png']),'-painters','-r360');
+end
+
+[sigu(1) mixx1(1) sigv(1) miyy1(1) sigw(1) mizz1(1)]
+%%
+% save([dirout,'histdata_withselfcal.mat'],'vecx1','vecx','vecy1','vecy','vecz1','vecz',...
+%     'Nex','Ney','Nez','Nimx','Nimy','Nimz','sigu','mixx1',...
+%     'sigv','miyy1','sigw','mizz1');
+
+close all;
+
+% % 
+% % vecx=linspace(0,ulx,Nb);
+% % vecy=linspace(0,uly,Nb);
+% % vecz=linspace(0,ulz,Nb);
+% % 
+% % Nex=histc(errx,vecx);
+% % Nimx=histc(unim_u,vecx);
+% % 
+% % Ney=histc(erry,vecy);
+% % Nimy=histc(unim_v,vecy);
+% % 
+% % Nez=histc(errz,vecz);
+% % Nimz=histc(unim_w,vecz);
+% % 
+% % 
+% % % X uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecx,Nex,'-k',vecx,Nimx,'r-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigu=rms(errx(errx<=ulx)).*ones(length(ll),1);
+% % mixx1=rms(unim_u(unim_u<=ulx)).*ones(length(ll),1);
+% % plot(sigu,ll,'k--',mixx1,ll,'r--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{x}',['\sigma_{x}^{',meth,'}']);axis([0 ulx 0 inf]);
+% % % title(['U Uncertainty Histogram coverage= ',num2str(covx1,'%02.02f')]);
+% % title('U Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Uhist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Uhist.png']),'-painters','-r360');
+% % end
+% % 
+% % % Y uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecy,Ney,'-k',vecy,Nimy,'b-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigv=rms(erry(erry<=uly)).*ones(length(ll),1);
+% % miyy1=rms(unim_v(unim_v<=uly)).*ones(length(ll),1);
+% % plot(sigv,ll,'k--',miyy1,ll,'b--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{y}',['\sigma_{y}^{',meth,'}']);axis([0 uly 0 inf]);
+% % % title(['V uncertainty histogram coverage= ',num2str(covy1,'%02.02f')]);
+% % title('V Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Vhist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Vhist.png']),'-painters','-r360');
+% % end
+% % 
+% % % Z uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecz,Nez,'-k',vecz,Nimz,'g-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigw=rms(errz(errz<=ulz)).*ones(length(ll),1);
+% % mizz1=rms(unim_w(unim_w<=ulz)).*ones(length(ll),1);
+% % plot(sigw,ll,'k--',mizz1,ll,'g--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{z}',['\sigma_{z}^{',meth,'}']);axis([0 ulz 0 inf]);
+% % % title(['W uncertainty histogram coverage= ',num2str(covz1,'%02.02f')]);
+% % title('W Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Whist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Whist.png']),'-painters','-r360');
+% % end
+% % 
+% % [sigu(1) mixx1(1) sigv(1) miyy1(1) sigw(1) mizz1(1)]
\ No newline at end of file
diff --git a/vortex_ring_postprocess_codes/Postprocess_VR_prana_01_11_2016.m b/vortex_ring_postprocess_codes/Postprocess_VR_prana_01_11_2016.m
new file mode 100644
index 0000000..949ca72
--- /dev/null
+++ b/vortex_ring_postprocess_codes/Postprocess_VR_prana_01_11_2016.m
@@ -0,0 +1,680 @@
+%Postprocess_stereo_uncertainty_vortexring thicksheet
+
+%Uncertainty Propagation for Stereo Geometric Reconstruction(Willert
+%Method)for vortex ring
+
+clear all;
+lw=4;fs=18;
+if ispc
+    main_dir='Z:\Stereo_Uncertainty_Work\';
+    addpath Z:\Stereo_Uncertainty_Work\Codes\New_codes\;
+
+else
+    main_dir='/home/shannon/a/bhattac3/Stereo_Uncertainty_Work/';
+    addpath /home/shannon/a/bhattac3/Stereo_Uncertainty_Work/Codes/New_codes/;
+
+end
+
+meth='IM';
+
+printfig=0;
+ensdispfldflag=0;
+
+wres=32;
+
+cam1='1';cam2='3';
+% cam1='2';cam2='4';
+
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    casename='IMcam13';
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    casename='IMcam24';
+end
+
+
+dirout=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','plots_final','Wihselfcal_angle');
+if ~exist(dirout,'dir')
+    mkdir(dirout);
+end
+dirout=[dirout,filesep,casename,'_'];
+
+
+stdjob= load(fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','jobs',['cam',cam1,cam2,'_VR_prana_fulljob_withselfcal.mat']));
+caljob=stdjob.caljobfile.datasave.caljob;
+calmat=[caljob.aXcam1 caljob.aYcam1 caljob.aXcam2 caljob.aYcam2];
+calmat;
+planarjob=stdjob.pranajob.Data;
+
+dispfield=load(fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Willert','dispfield_ensemble',['VR',cam1,cam2,'_scc_pass2_00025.mat']));
+
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    imagedir.imdir1=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images1','E_camera_1_frame_');
+    imagedir.imdir2=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images2','E_camera_3_frame_');
+
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    imagedir.imdir1=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images1','Camera_02_frame_');
+    imagedir.imdir2=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Dewarped Images2','Camera_04_frame_');
+
+end
+
+%% Input Directories
+basedir=fullfile(main_dir,'New_Stereo_uncertainty_tests','Vortex_Ring','Prana_process','results',['cam',cam1,cam2],'With_selfcal','Willert',filesep);
+
+dirwc1=[basedir,'dispfields',filesep,'VR',cam1,cam2,'_scc_pass4_'];
+istring1=sprintf(['%%s%%0%0.0fd.','mat'],5);
+istring2=sprintf(['%%s%%0%0.0fd.','mat'],5);
+
+%Camera1 directory
+dircam1a=[basedir,'Camera',cam1,filesep,'VR',cam1,cam2,'_scc_pass4_'];
+
+%Camera2 directory
+dircam2a=[basedir,'Camera',cam2,filesep,'VR',cam1,cam2,'_scc_pass4_'];
+
+%reconstructed field directories
+dirrec=[basedir,'Camera',cam1,'Camera',cam2,'_3Cfields',filesep,'piv_2d3c_cam',cam1,'cam',cam2,'_pass_4_'];
+
+%true solution directory
+dirtrue=fullfile(main_dir,'VortexRing_True_solution','MATLAB_Data','Eval_1','E_eval1_TrueSolution_');
+%dirtrue='/home/shannon/a/bhattac3/Stereo_Uncertainty_Work/Process_vortexring_cam2cam4/Velocity_field_with_selfcal/Willert/Camera2Camera4_3Cfields/piv_2d3c_cam2cam4_pass_4_';
+
+scaling=stdjob.scaling.wil;
+mx=scaling.xscale;
+my=scaling.yscale;
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+    mz=scaling.yscale;
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+    mz=scaling.xscale;
+end
+%seconds per frame
+spf=0.001;
+
+t=50;   
+%get uncertainties in x and y
+ofs=24;
+%% Find uncertainty in angles  
+outputdir=dirout;
+fstart=ofs;
+Nim=t;
+
+[Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,alpha1,alpha2,beta1,beta2,Umx,Umy]=stereo_calcoeff_uncertainty(caljob,planarjob,dispfield,imagedir,outputdir,ensdispfldflag,wres,fstart,Nim);
+ 
+% [mean(alpha1(:)) mean(alpha2(:)) mean(beta1(:)) mean(beta2(:))]
+
+% [mean(rad2deg(alpha1(:))) mean(rad2deg(alpha2(:))) mean(rad2deg(beta1(:))) mean(rad2deg(beta2(:)))]
+% keyboard; 
+% [mean(Un_alpha1(:)) mean(Un_alpha2(:)) mean(Un_beta1(:)) mean(Un_beta2(:))]
+
+%%  
+[Sx,Sy]=size(alpha1);
+
+for i=1:t
+  
+   %% Load reconstructed solution
+   wilsol=load(sprintf(istring2,dirrec,i+ofs)); 
+   Us(:,:,i)=(spf/mx)*1e3*wilsol.U;
+   Vs(:,:,i)=(spf/my)*1e3*wilsol.V;
+   Ws(:,:,i)=(spf/mz)*1e3*wilsol.W;
+
+   xgrid=1000*wilsol.X;
+   ygrid=1000*wilsol.Y;
+   zgrid=1000*wilsol.Z;
+   
+   [r,c]=size(xgrid);
+   
+   %%  Load true solution
+   TS1=load(sprintf(istring2,dirtrue,i+24));
+   TS=TS1.velocity_data;
+   ut=(spf/mx)*1000*TS.x_velocity;
+   vt=(spf/my)*1000*TS.y_velocity;
+   wt=(spf/mz)*1000*TS.z_velocity;
+   [xt,yt]=meshgrid(-34:0.45:20,-27:0.45:27);
+
+%    Ut(:,:,i)=interp2(xt,yt,ut,xgrid,ygrid,'spline',0);
+%    Vt(:,:,i)=interp2(xt,yt,vt,xgrid,ygrid,'spline',0);
+%    Wt(:,:,i)=interp2(xt,yt,wt,xgrid,ygrid,'spline',0);
+   Ut(:,:,i)=interp2(xt,yt,ut,xgrid,ygrid,'linear',0);
+   Vt(:,:,i)=interp2(xt,yt,vt,xgrid,ygrid,'linear',0);
+   Wt(:,:,i)=interp2(xt,yt,wt,xgrid,ygrid,'linear',0);
+
+% keyboard;
+% figure;pcolor(xt,yt,wt);shading interp;hold on;
+% axis image;
+% axis off;
+% skp=4;
+% quiver(xt(1:skp:end,1:skp:end),yt(1:skp:end,1:skp:end),ut(1:skp:end,1:skp:end),vt(1:skp:end,1:skp:end),1,'black');
+% print(gcf,'-dpng','vortex_ring2_frame2.png','-r300');
+   
+%    figure(4);
+%    subplot(2,3,1);imagesc(xgrid(:),ygrid(:),Us(:,:,i));colorbar;title('Us');axis equal square;caxis([-1 3]);
+%    subplot(2,3,2);imagesc(xgrid(:),ygrid(:),Vs(:,:,i));colorbar;title('Vs');axis equal square;caxis([-1.5 1.5]);
+%    subplot(2,3,3);imagesc(xgrid(:),ygrid(:),Ws(:,:,i));colorbar;title('Ws');axis equal square;caxis([-1 1]);
+%    subplot(2,3,4);imagesc(xgrid(:),ygrid(:),Ut(:,:,i));colorbar;title('Ut');axis equal square;caxis([-1 3]);
+%    subplot(2,3,5);imagesc(xgrid(:),ygrid(:),Vt(:,:,i));colorbar;title('Vt');axis equal square;caxis([-1.5 1.5]);
+%    subplot(2,3,6);imagesc(xgrid(:),ygrid(:),Wt(:,:,i));colorbar;title('Wt');axis equal square;caxis([-1 1]);
+%    keyboard;
+% %    pause(0.05);
+   
+   %% Get error fields
+   errx1(:,:,i)=(Ut(:,:,i)-Us(:,:,i));
+   erry1(:,:,i)=(Vt(:,:,i)-Vs(:,:,i));
+   errz1(:,:,i)=(Wt(:,:,i)-Ws(:,:,i));
+   
+   uttemp=Ut(:,:,i);
+   vttemp=Vt(:,:,i);
+   ustemp=Us(:,:,i);
+   vstemp=Vs(:,:,i);
+   
+   truediv(:,:,i)=divergence(xgrid,ygrid,uttemp,vttemp);
+   soldiv(:,:,i)=divergence(xgrid,ygrid,ustemp,vstemp);
+   
+   %% Load cam1cam2 2d fields
+   cam1sola=load(sprintf(istring2,dircam1a,i+ofs)); 
+   cam2sola=load(sprintf(istring2,dircam2a,i+ofs)); 
+    
+   %get 2d velocities
+   U1=cam1sola.U;
+   V1=cam1sola.V;
+   X1=cam1sola.X;
+   Y1=cam1sola.Y;
+   
+   U2=cam2sola.U;
+   V2=cam2sola.V;
+   X2=cam2sola.X;
+   Y2=cam2sola.Y;
+   
+   %% Get 2d uncertainties
+   %Camera1
+   imx1(:,:,i)=cam1sola.imx;
+   imy1(:,:,i)=cam1sola.imy;
+   Nim1(:,:,i)=cam1sola.Npart;
+   
+   %Camera2
+   imx2(:,:,i)=cam2sola.imx;
+   imy2(:,:,i)=cam2sola.imy;
+   Nim2(:,:,i)=cam2sola.Npart;
+  
+  
+   %% Stereo Uncertainty Propagation
+
+    Unu1=imx1(:,:,i);
+    Unu2=imx2(:,:,i);
+    Unv1=imy1(:,:,i);
+    Unv2=imy2(:,:,i);
+
+   [Un_u1,Un_v1,Un_w1,JU,JV,JW]=stereo_uncertainty_propagation_withcalcoeff(Unu1,Unv1,Unu2,Unv2,U1,V1,U2,V2,Ws(:,:,i),Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,alpha1,alpha2,beta1,beta2,mx,my);
+
+   [Un_u3,Un_v3,Un_w3]=stereo_uncertainty_propagation_noangle(Unu1,Unv1,Unu2,Unv2,alpha1,alpha2,beta1,beta2,mx,my);
+   
+   % Uncertainty in magnification
+%    Umagx(:,:,i)=sqrt((Umx^2)*((Us(:,:,i)*mx).^2)/(mx^4));
+%    Umagy(:,:,i)=sqrt((Umy^2)*((Vs(:,:,i)*my).^2)/(my^4));
+%    Umagz(:,:,i)=sqrt((Umx^2)*((Ws(:,:,i)*mx).^2)/(mx^4));
+   
+%     unim_u(:,:,i)=sqrt(Un_u1.^2+Umagx(:,:,i).^2);
+%     unim_v(:,:,i)=sqrt(Un_v1.^2+Umagy(:,:,i).^2);
+%     unim_w(:,:,i)=sqrt(Un_w1.^2+Umagz(:,:,i).^2);
+%     unim_um(:,:,i)=Un_u3;
+%     unim_vm(:,:,i)=Un_v3;
+%     unim_wm(:,:,i)=Un_w3;
+   
+    unim_u(:,:,i)=Un_u1;
+    unim_v(:,:,i)=Un_v1;
+    unim_w(:,:,i)=Un_w1;
+    unim_um(:,:,i)=Un_u3;
+    unim_vm(:,:,i)=Un_v3;
+    unim_wm(:,:,i)=Un_w3;
+
+%     unim_u(:,:,i)=Un_u3;
+%     unim_v(:,:,i)=Un_v3;
+%     unim_w(:,:,i)=Un_w3;
+%     unim_um(:,:,i)=Un_u1;
+%     unim_vm(:,:,i)=Un_v1;
+%     unim_wm(:,:,i)=Un_w1;
+
+    
+
+end
+soldiv2=soldiv(abs(soldiv)<1);
+[mean(truediv(:)) mean(soldiv2(:))]
+[rms(truediv(:)) rms(soldiv2(:))]
+keyboard;
+%Good measurement cut off
+ulx=0.5;uly=0.5;ulz=1.5;
+
+%% Plot RMS Uncertainty and Error fields (Spatial RMS Profile)
+
+for i=1:Sx
+    for j=1:Sy
+        
+        
+        rms2derrx=squeeze(abs(errx1(i,j,:)));
+        rms2dunx=squeeze(abs(unim_u(i,j,:)));
+        ix1=find(rms2derrx>ulx);
+        rms2derrx(ix1)=[];
+        rms2dunx(ix1)=[];
+%         Nx(i,j)=length(ix1);
+        
+        rms2derry=squeeze(abs(erry1(i,j,:)));
+        rms2duny=squeeze(abs(unim_v(i,j,:)));
+        iy1=find(rms2derry>uly);
+        rms2derry(iy1)=[];
+        rms2duny(iy1)=[];
+%         Ny(i,j)=length(iy1);
+        
+        rms2derrz=squeeze(abs(errz1(i,j,:)));
+        rms2dunz=squeeze(abs(unim_w(i,j,:)));
+        iz1=find(rms2derrz>ulz);
+        rms2derrz(iz1)=[];
+        rms2dunz(iz1)=[];
+%         Nz(i,j)=length(iz1);
+       
+        
+        Ue(i,j)=rms(rms2derrx);
+        Uun(i,j)=rms(rms2dunx);
+        Ve(i,j)=rms(rms2derry);
+        Vun(i,j)=rms(rms2duny);
+        We(i,j)=rms(rms2derrz);
+        Wun(i,j)=rms(rms2dunz);
+        
+    end
+end
+figure;
+subplot(2,3,1);imagesc(Ue);caxis([0 0.2]);
+subplot(2,3,2);imagesc(Ve);caxis([0 0.2]);
+subplot(2,3,3);imagesc(We);caxis([0 0.4]);
+subplot(2,3,4);imagesc(Uun);caxis([0 0.2]);
+subplot(2,3,5);imagesc(Vun);caxis([0 0.2]);
+subplot(2,3,6);imagesc(Wun);caxis([0 0.4]);
+
+%% Plot RMS Uncertainty and Error fields (Temporal RMS Profile)
+
+for k=1:t
+    
+        rex=squeeze(abs(errx1(:,:,k)));
+        rIMx=squeeze(abs(unim_u(:,:,k)));
+        rex=rex(:);
+        rIMx=rIMx(:);
+        ix1=find(rex>ulx);
+        rex(ix1)=[];
+        rIMx(ix1)=[];
+%         Nx(i,j)=length(ix1);
+        
+        rey=squeeze(abs(erry1(:,:,k)));
+        rIMy=squeeze(abs(unim_v(:,:,k)));
+        rey=rey(:);
+        rIMy=rIMy(:);
+        iy1=find(rey>uly);
+        rey(iy1)=[];
+        rIMy(iy1)=[];
+%         Ny(i,j)=length(iy1);
+        
+        rez=squeeze(abs(errz1(:,:,k)));
+        rIMz=squeeze(abs(unim_w(:,:,k)));
+        rez=rez(:);
+        rIMz=rIMz(:);
+        iz1=find(rez>ulz);
+        rez(iz1)=[];
+        rIMz(iz1)=[];
+%         Nz(i,j)=length(iz1);
+       
+        
+        siguerr(k)=rms(rex);
+        siguun(k)=rms(rIMx);
+        sigverr(k)=rms(rey);
+        sigvun(k)=rms(rIMy);
+        sigwerr(k)=rms(rez);
+        sigwun(k)=rms(rIMz);
+        
+   
+end
+
+figure;hold on;
+plot(1:t,siguerr,'r-o',1:t,siguun,'r-*');
+plot(1:t,sigverr,'b-o',1:t,sigvun,'b-*');
+plot(1:t,sigwerr,'g-o',1:t,sigwun,'g-*');
+hold off;
+title('RMS Time Series (Cam24)');xlabel('Frame No.');
+ylabel('RMS error, uncertainty(pix/frame)');
+% hleg=legend({'\sigma^{e}_{u}','\sigma^{CS}_{u}','\sigma^{e}_{v}','\sigma^{CS}_{v}','\sigma^{e}_{w}','\sigma^{CS}_{w}'});
+% set(hleg,'location','eastoutside');
+set(gca,'FontSize',20);
+axis([0 50 0 0.5]);axis square;
+%% RMS 2D Uncertainty spatial maps
+nc=20;
+IMX1=rms(imx1,3);
+IMY1=rms(imy1,3);
+IMX2=rms(imx2,3);
+IMY2=rms(imy2,3);
+figure;
+subplot(2,2,1);contour(xgrid,ygrid,IMX1,nc);caxis([0 0.3]);title('\sigma^{IM}_{U1}');axis square;
+subplot(2,2,2);contour(xgrid,ygrid,IMY1,nc);caxis([0 0.3]);title('\sigma^{IM}_{V1}');axis square;
+subplot(2,2,3);contour(xgrid,ygrid,IMX2,nc);caxis([0 0.3]);title('\sigma^{IM}_{U2}');axis square;
+subplot(2,2,4);contour(xgrid,ygrid,IMY2,nc);caxis([0 0.3]);title('\sigma^{IM}_{V2}');axis square;
+
+%%
+save([meth,'_Cam',cam1,cam2,'_RMSprofiles.mat'],'xgrid','ygrid','IMX1',...
+    'IMX2','IMY1','IMY2','siguerr','siguun','sigverr','sigvun','sigwerr','sigwun',...
+    'Ue','Uun','Ve','Vun','We','Wun');
+keyboard;
+%%
+%[mean(errx(:)) mean(erry(:)) mean(errz(:))]
+
+% [mean(imx1(:)) mean(imy1(:)) mean(Nim1(:)) mean(imx2(:)) mean(imy2(:)) mean(Nim2(:))]
+% [rms(imx1(:)) rms(imy1(:)) mean(Nim1(:)) rms(imx2(:)) rms(imy2(:)) mean(Nim2(:))]
+
+%%
+errx=(errx1);
+erry=(erry1);
+errz=(errz1);
+
+% errx1(abs(errx1)>0.3)=[];
+% erry1(abs(erry1)>0.3)=[];
+% errz1(abs(errz1)>0.3)=[];
+% 
+% Nbin1=60;
+% figure;hist(errx1(:),Nbin1);
+% figure;hist(erry1(:),Nbin1);
+% figure;hist(errz1(:),Nbin1);
+
+Nc=100;
+
+ul=0.5;
+
+vx1=linspace(0,ul,Nc);
+vy1=linspace(0,ul,Nc);
+
+vx2=linspace(0,ul,Nc);
+vy2=linspace(0,ul,Nc);
+
+Nu1=histc(imx1(:),vx1);
+Nv1=histc(imy1(:),vy1);
+Nu2=histc(imx2(:),vx2);
+Nv2=histc(imy2(:),vy2);
+
+figure;set(gcf,'DefaultLineLineWidth',1);set(gca,'FontSize',fs);
+plot(vx1,Nu1,'r-o',vy1,Nv1,'b-o',vx2,Nu2,'r-+',vy2,Nv2,'b-+');
+title('Planar Correlation Uncertainties');
+legend('U1','V1','U2','V2');
+if printfig==1
+    print(gcf,'-dpng',[dirout,'planarhist.png'],'-r300');
+end
+%% RMS error and uncertainty field
+
+% rex=rms(errx,3);
+% rey=rms(erry,3);
+% rez=rms(errz,3);
+% rIMx=rms(unim_u,3);
+% rIMy=rms(unim_v,3);
+% rIMz=rms(unim_w,3);
+% 
+% figure;imagesc(rex(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{e}_{u}');axis square tight;
+% figure;imagesc(rIMx(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{CS}_{u}');axis square tight;
+% figure;imagesc(rez(end:-1:1,:));caxis([0 0.5]);axis off;title('\sigma^{e}_{w}');axis square tight;
+% figure;imagesc(rIMz(end:-1:1,:));caxis([0 0.5]);axis off;title('\sigma^{CS}_{w}');axis square tight;
+% figure;imagesc(rey(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{e}_{v}');axis square tight;
+% figure;imagesc(rIMy(end:-1:1,:));caxis([0 0.25]);axis off;title('\sigma^{CS}_{v}');axis square tight;
+% keyboard;
+%% eliminate invalid measurements
+errx=errx(:);
+erry=erry(:);
+errz=errz(:);
+unim_u=unim_u(:);
+unim_v=unim_v(:);
+unim_w=unim_w(:);
+
+
+unim_um=unim_um(:);
+unim_vm=unim_vm(:);
+unim_wm=unim_wm(:);
+
+% %Good measurement cut off
+% ulx=0.5;uly=0.5;ulz=1.5;
+
+inx=find(abs(errx(:))>ulx);
+iny=find(abs(erry(:))>uly);
+inz=find(abs(errz(:))>ulz);
+
+errx(inx)=[];
+erry(iny)=[];
+errz(inz)=[];
+unim_u(inx)=[];
+unim_v(iny)=[];
+unim_w(inz)=[];
+
+
+unim_um(inx)=[];
+unim_vm(iny)=[];
+unim_wm(inz)=[];
+
+% [mean(errx(:)) mean(erry(:)) mean(errz(:))]
+%%
+
+% Calculate coverage
+cnt1=0;cnt2=0;cnt3=0;cnt4=0;cnt5=0;cnt6=0;
+l1=length(errx);
+l2=length(erry);
+l3=length(errz);
+for m1=1:l1
+    if abs(errx(m1))<=unim_u(m1)
+        cnt1=cnt1+1;
+    end
+    if abs(errx(m1))<=unim_um(m1)
+        cnt4=cnt4+1;
+    end
+    
+end
+
+for m2=1:l2
+    if abs(erry(m2))<=unim_v(m2)
+        cnt2=cnt2+1;
+    end
+    if abs(erry(m2))<=unim_vm(m2)
+        cnt5=cnt5+1;
+    end
+end
+
+for m3=1:l3
+    
+    if abs(errz(m3))<=unim_w(m3)
+        cnt3=cnt3+1;
+    end
+    if abs(errz(m3))<=unim_wm(m3)
+        cnt6=cnt6+1;
+    end
+    
+end
+covx1=100*cnt1/l1;
+covy1=100*cnt2/l2;
+covz1=100*cnt3/l3;
+fprintf('\n Coverage with angle uncertainty');
+[covx1 covy1 covz1]
+covx2=100*cnt4/l1;
+covy2=100*cnt5/l2;
+covz2=100*cnt6/l3;
+fprintf('\n Coverage without angle uncertainty');
+[covx2 covy2 covz2]
+
+%%
+%% Plot Histograms
+if ispc
+addpath Z:\Planar_Uncertainty_work\codes\export_fig\;
+addpath Z:\Planar_Uncertainty_work\codes\color_code\;
+else
+    addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/export_fig/;
+    addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/color_code/;
+end
+
+set(0,'DefaultAxesFontName', 'Times New Roman');
+
+%Plot histograms
+Nb=30;
+
+vecx=linspace(0,ulx,Nb);
+vecy=linspace(0,uly,Nb);
+vecz=linspace(0,ulz,Nb);
+vecx1=linspace(-ulx,ulx,2*Nb-1);
+vecy1=linspace(-uly,uly,2*Nb-1);
+vecz1=linspace(-ulz,ulz,2*Nb-1);
+
+Nex=histc(errx,vecx1)/length(errx);
+Nimx=histc(unim_u,vecx)/length(errx);
+
+Ney=histc(erry,vecy1)/length(erry);
+Nimy=histc(unim_v,vecy)/length(erry);
+
+Nez=histc(errz,vecz1)/length(errz);
+Nimz=histc(unim_w,vecz)/length(errz);
+
+
+% X uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecx1,Nex,'-k');
+plot(vecx,Nimx,'r-');
+% HL1=fill(vecx1,Nex,'-k');
+% HL2=fill(vecx,Nimx,'r-');
+% set(HL1,'facealpha',.5);
+% set(HL2,'facealpha',.5);
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigu=rms(errx(abs(errx)<=ulx)).*ones(length(ll),1);
+mixx1=rms(unim_u(unim_u<=ulx)).*ones(length(ll),1);
+plot(sigu,ll,'k--',mixx1,ll,'r--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{x}',['\sigma_{x}^{',meth,'}']);
+axis([-ulx ulx 0 inf]);
+% title(['U Uncertainty Histogram coverage= ',num2str(covx1,'%02.02f')]);
+title('U Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Uhist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Uhist.png']),'-painters','-r360');
+end
+
+% Y uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecy1,Ney,'-k');
+plot(vecy,Nimy,'b-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigv=rms(erry(abs(erry)<=uly)).*ones(length(ll),1);
+miyy1=rms(unim_v(unim_v<=uly)).*ones(length(ll),1);
+plot(sigv,ll,'k--',miyy1,ll,'b--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{y}',['\sigma_{y}^{',meth,'}']);
+axis([-uly uly 0 inf]);
+% title(['V uncertainty histogram coverage= ',num2str(covy1,'%02.02f')]);
+title('V Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Vhist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Vhist.png']),'-painters','-r360');
+end
+
+% Z uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecz1,Nez,'-k');
+plot(vecz,Nimz,'g-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigw=rms(errz(abs(errz)<=ulz)).*ones(length(ll),1);
+mizz1=rms(unim_w(unim_w<=ulz)).*ones(length(ll),1);
+plot(sigw,ll,'k--',mizz1,ll,'g--'); 
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{z}',['\sigma_{z}^{',meth,'}']);
+axis([-ulz ulz 0 inf]);
+% title(['W uncertainty histogram coverage= ',num2str(covz1,'%02.02f')]);
+title('W Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+%     print(gcf,'-dpng',[dirout,'Whist.png'],'-r300');
+    export_fig(gcf,fullfile([dirout,'Whist.png']),'-painters','-r360');
+end
+
+[sigu(1) mixx1(1) sigv(1) miyy1(1) sigw(1) mizz1(1)]
+%%
+% save([dirout,'histdata_withselfcal.mat'],'vecx1','vecx','vecy1','vecy','vecz1','vecz',...
+%     'Nex','Ney','Nez','Nimx','Nimy','Nimz','sigu','mixx1',...
+%     'sigv','miyy1','sigw','mizz1');
+% 
+% close all;
+
+%%
+
+% % vecx=linspace(0,ulx,Nb);
+% % vecy=linspace(0,uly,Nb);
+% % vecz=linspace(0,ulz,Nb);
+% % 
+% % Nex=histc(errx,vecx);
+% % Nimx=histc(unim_u,vecx);
+% % 
+% % Ney=histc(erry,vecy);
+% % Nimy=histc(unim_v,vecy);
+% % 
+% % Nez=histc(errz,vecz);
+% % Nimz=histc(unim_w,vecz);
+% % 
+% % 
+% % % X uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecx,Nex,'-k',vecx,Nimx,'r-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigu=rms(errx(errx<=ulx)).*ones(length(ll),1);
+% % mixx1=rms(unim_u(unim_u<=ulx)).*ones(length(ll),1);
+% % plot(sigu,ll,'k--',mixx1,ll,'r--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{x}',['\sigma_{x}^{',meth,'}']);axis([0 ulx 0 inf]);
+% % % title(['U Uncertainty Histogram coverage= ',num2str(covx1,'%02.02f')]);
+% % title('U Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Uhist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Uhist.png']),'-painters','-r360');
+% % end
+% % 
+% % % Y uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecy,Ney,'-k',vecy,Nimy,'b-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigv=rms(erry(erry<=uly)).*ones(length(ll),1);
+% % miyy1=rms(unim_v(unim_v<=uly)).*ones(length(ll),1);
+% % plot(sigv,ll,'k--',miyy1,ll,'b--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{y}',['\sigma_{y}^{',meth,'}']);axis([0 uly 0 inf]);
+% % % title(['V uncertainty histogram coverage= ',num2str(covy1,'%02.02f')]);
+% % title('V Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Vhist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Vhist.png']),'-painters','-r360');
+% % end
+% % 
+% % % Z uncertainty Histogram
+% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+% % plot(vecz,Nez,'-k',vecz,Nimz,'g-');
+% % yl=ylim(gca);ll=linspace(1,max(yl(2)),10);
+% % sigw=rms(errz(errz<=ulz)).*ones(length(ll),1);
+% % mizz1=rms(unim_w(unim_w<=ulz)).*ones(length(ll),1);
+% % plot(sigw,ll,'k--',mizz1,ll,'g--'); 
+% % 
+% % grid on;box on; set(gcf,'color','white');hold off;
+% % 
+% % legend('e_{z}',['\sigma_{z}^{',meth,'}']);axis([0 ulz 0 inf]);
+% % % title(['W uncertainty histogram coverage= ',num2str(covz1,'%02.02f')]);
+% % title('W Uncertainty Histogram');
+% % xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+% % if printfig==1
+% % %     print(gcf,'-dpng',[dirout,'Whist.png'],'-r300');
+% %     export_fig(gcf,fullfile([dirout,'Whist.png']),'-painters','-r360');
+% % end
+% % 
+% % [sigu(1) mixx1(1) sigv(1) miyy1(1) sigw(1) mizz1(1)]
\ No newline at end of file