%Compare motion correction using slightly different algorithms (mailny
%filter cutoff frequencies and integration techniques) between PSD and
%UConn
%Selected a calm day (doy320) and the stormy day with more waves (doy328)  

clear all;fclose all;
rdcon=180/pi;
skipit=0; %to skip direct comparison plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%selects the day to compare
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ddd=320; %bad reldir that day! Comparison not that great!
% ddd=316; 
ddd=328; 
% ddd=335; 
jd=sprintf('%03i',ddd);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%brings in 20-min chunks UConn data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ydayS: Day-Of-Year (GPS sync.)
% euler: euler angles (UConn system)
% eulerE: euler angles (PSD system)
% uvwplat: platform motion with ship's speed (UConn system)
% xyzplat: platform displacements with ship's speed (UConn system)
% uvwplatE: platform motion with ship's speed (PSD system)
% xyzplatE: platform displacements with ship's speed (PSD system)
% uvwplatmot: platform motion only (UConn system)
% uvwplatmotE: platform motion only (PSD system)
% uvwmot: corrected sonic velocities without GPS (UConn system)
% uvwmotE: corrected sonic velocities without GPS (PSD system)
% uvw: corrected sonic velocities with GPS (UConn system) --> was done 
% using high resolution GPS rather than mean
% uvwE: corrected sonic velocities with GPS (PSD system)
% hr20=[hhh*3600 (hhh+20/60)*3600 (hhh+40/60)*3600];
% jdUConn=[];UConn_euler=[];UConn_eulerPSD=[];UConn_uvwplatmot=[];UConn_uvwplatmotPSD=[];
% for ii=1:length(hr20)
%     UConnmot=load(['E:\DYNAMO_2011\Revelle\flux\Raw\UConn data\PSDUConn_motcorr_comparison\2011' jd sprintf('%05g',hr20(ii)) '.mat']);
%     jdUConn=[jdUConn UConnmot.ydayS];
%     UConn_euler=[UConn_euler UConnmot.euler];
%     UConn_eulerPSD=[UConn_eulerPSD UConnmot.eulerE];
%     UConn_uvwplatmot=[UConn_uvwplatmot UConnmot.uvwplatmot];
%     UConn_uvwplatmotPSD=[UConn_uvwplatmotPSD UConnmot.uvwplatmotE];
% end
% clear UConnmot
%%-->instead of concatenating one hour, Jim suggests to compare 20min
%%chunks as to not introduce a 1 pt lag here and there (distortion)
dr=dir(['E:\DYNAMO_2011\Revelle\flux\Raw\UConn data\PSDUConn_motcorr_comparison\2011' jd '*.mat']);
specview=ones(length(dr),54*12);

for ii=1:length(dr)
% for ii=38; %1:length(dr)
    filestub=dr(ii).name(1:12);
    UConnmot=load(['E:\DYNAMO_2011\Revelle\flux\Raw\UConn data\PSDUConn_motcorr_comparison\' filestub '.mat']);
    ti=str2num(filestub(end-4:end)); 
    jdUConn=UConnmot.ydayS;
    UConn_euler=UConnmot.euler*rdcon;
    UConn_eulerPSD=UConnmot.eulerE*rdcon;
    UConn_uvwplatmot=UConnmot.uvwplatmot;
    UConn_uvwplatmotPSD=UConnmot.uvwplatmotE;
    UConn_uvwmot=UConnmot.uvwmot;
    UConn_uvwmotPSD=UConnmot.uvwmotE;
    UConn_uvw=UConnmot.uvw;
    UConn_uvwPSD=UConnmot.uvwE;
% end
[Uj,Vj,Wj,alphaj,betaj]=rotcorr(UConn_uvw(1,:),UConn_uvw(2,:),UConn_uvw(3,:));    
clear UConn_uvw   
UConn_uvw=[Uj;Vj;Wj];
clear temp Uj Vj Wj alphaj betaj
[Uj,Vj,Wj,alphaj,betaj]=rotcorr(UConn_uvwPSD(1,:),UConn_uvwPSD(2,:),UConn_uvwPSD(3,:));    
clear UConn_uvwPSD   
UConn_uvwPSD=[Uj;Vj;Wj];
clear temp Uj Vj Wj alphaj betaj
[Uj,Vj,Wj,alphaj,betaj]=rotcorr(UConn_uvwmot(1,:),UConn_uvwmot(2,:),UConn_uvwmot(3,:));    
clear UConn_uvwmot  
UConn_uvwmot=[Uj;Vj;Wj];
clear temp Uj Vj Wj alphaj betaj
[Uj,Vj,Wj,alphaj,betaj]=rotcorr(UConn_uvwmotPSD(1,:),UConn_uvwmotPSD(2,:),UConn_uvwmotPSD(3,:));    
clear UConn_uvwmotPSD  
UConn_uvwmotPSD=[Uj;Vj;Wj];
clear temp Uj Vj Wj alphaj betaj

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%brings in PSD hourly data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hhh=floor(ti/3600);hr=sprintf('%02i',hhh);  
PSDmot=load(['E:\Desktop_bd_Cdrive\Data\CRUISES\DYNAMO_2011\flux_data\MotionData\DYNAMO2011_motionfile_',jd hr,'.prn']);
% motmatall=[tref' uvwplat(1,:)' uvwplat(2,:)' uvwplat(3,:)' xyzplat(1,:)' xyzplat(2,:)' xyzplat(3,:)' euler(1,:)' euler(2,:)' euler(3,:)' deg_rate(1,:)' deg_rate(2,:)' deg_rate(3,:)'];
jdPSD=ddd+hhh/24+PSDmot(:,1)/86400;
llPSD=find(jdPSD>=jdUConn(1) & jdPSD<jdUConn(end));
PSD_uvwplatmot=[PSDmot(llPSD,2) PSDmot(llPSD,3) PSDmot(llPSD,4)];
PSD_xyzplatmot=[PSDmot(llPSD,5) PSDmot(llPSD,6) PSDmot(llPSD,7)];
PSD_euler=[PSDmot(llPSD,8) PSDmot(llPSD,9) PSDmot(llPSD,10)]*rdcon;

PSDuvwmat=load(['E:\Desktop_bd_Cdrive\Data\CRUISES\DYNAMO_2011\flux_data\MotionData\DYNAMO2011_UVWstream_',jd hr,'.prn ']);
PSD_uvw=[PSDuvwmat(llPSD,2) PSDuvwmat(llPSD,3) PSDuvwmat(llPSD,4)];

if skipit
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %interpolates
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    PSD_uvwplatmoti=interp1(jdPSD(llPSD),PSD_uvwplatmot,jdUConn)';
    PSD_uvwi=interp1(jdPSD(llPSD),PSD_uvw,jdUConn)';
    PSD_euleri=interp1(jdPSD(llPSD),PSD_euler,jdUConn)';


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %finds possible lags and do comparison plots
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Compares corrected motion platform variables

    Xlag = findlag(UConn_uvwplatmot(3,:),PSD_uvwplatmoti(3,:),50,0);
    XlagPSD = findlag(UConn_uvwplatmotPSD(3,:),PSD_uvwplatmoti(3,:),50,0);
    figure(1);
    subplot(3,1,1);plot(jdUConn,PSD_uvwplatmoti(1,:),'b.--',jdUConn-abs(Xlag/10)/86400,UConn_uvwplatmot(1,:),'g*--',jdUConn-abs(XlagPSD/10)/86400,UConn_uvwplatmotPSD(1,:),'m*--')
    grid;xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('U platmot(m/s)');legend('PSD','UConn','UConnPSD')
    subplot(3,1,2);plot(jdUConn,PSD_uvwplatmoti(2,:),'b.--',jdUConn-abs(Xlag/10)/86400,UConn_uvwplatmot(2,:),'g*--',jdUConn-abs(XlagPSD/10)/86400,UConn_uvwplatmotPSD(2,:),'m*--')
    grid;xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('V platmot(m/s)');legend('PSD','UConn','UConnPSD')
    subplot(3,1,3);plot(jdUConn,PSD_uvwplatmoti(3,:),'b.--',jdUConn-abs(Xlag/10)/86400,UConn_uvwplatmot(3,:),'g*--',jdUConn-abs(XlagPSD/10)/86400,UConn_uvwplatmotPSD(3,:),'m*--')
    grid;xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('W platmot(m/s)');legend('PSD','UConn','UConnPSD')

    Xlagroll = findlag(UConn_euler(1,:),PSD_euleri(1,:),50,0)
    Xlagpitch = findlag(UConn_euler(2,:),PSD_euleri(2,:),50,0)
    Xlagyaw = findlag(UConn_euler(3,:),PSD_euleri(3,:),50,0)
    figure(2);
    subplot(3,1,1);plot(jdUConn,PSD_euleri(1,:),'b.--',jdUConn-abs(Xlagroll/10)/86400,UConn_euler(1,:),'g*--',jdUConn-abs(Xlagroll/10)/86400,UConn_eulerPSD(1,:),'m*--')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('roll (deg)');legend('PSD','UConn','UConnPSD')
    subplot(3,1,2);plot(jdUConn,PSD_euleri(2,:),'b.--',jdUConn-abs(Xlagpitch/10)/86400,UConn_euler(2,:),'g*--',jdUConn-abs(Xlagpitch/10)/86400,UConn_eulerPSD(2,:),'m*--')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('pitch (deg)');legend('PSD','UConn','UConnPSD')
    subplot(3,1,3);plot(jdUConn,PSD_euleri(3,:),'b.--',jdUConn-abs(Xlagyaw/10)/86400,UConn_euler(3,:),'g*--',jdUConn-abs(Xlagyaw/10)/86400,UConn_eulerPSD(3,:),'m*--')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('yaw (deg)');legend('PSD','UConn','UConnPSD')


    %removes the beginning and end of time series due to tapering filter
    %effects
    ptcut=200;
    figure(3);
    subplot(2,2,1);plot(PSD_uvwplatmoti(1,ptcut:end+Xlag),UConn_uvwplatmot(1,ptcut-Xlag:end),'.',PSD_uvwplatmoti(1,ptcut:end+XlagPSD),UConn_uvwplatmotPSD(1,ptcut-XlagPSD:end),'g.',[-0.4 0.4],[-0.4 0.4],'k--','LineWidth',2);
    grid;axis([-0.4 0.4 -0.4 0.4]);xlabel('PSD U platmot(m/s)');
    legend('UConn vs PSD','UConnPSD vs PSD','Location','NorthWest')
    subplot(2,2,2);plot(PSD_uvwplatmoti(2,ptcut:end+Xlag),UConn_uvwplatmot(2,ptcut-Xlag:end),'.',PSD_uvwplatmoti(2,ptcut:end+XlagPSD),UConn_uvwplatmotPSD(2,ptcut-XlagPSD:end),'g.',[-0.4 0.4],[-0.4 0.4],'k--','LineWidth',2);
    grid;axis([-0.4 0.4 -0.4 0.4]);xlabel('PSD V platmot(m/s)');
    legend('UConn vs PSD','UConnPSD vs PSD','Location','NorthWest')
    subplot(2,2,3);plot(PSD_uvwplatmoti(3,ptcut:end+Xlag),UConn_uvwplatmot(3,ptcut-Xlag:end),'.',PSD_uvwplatmoti(3,ptcut:end+XlagPSD),UConn_uvwplatmotPSD(3,ptcut-XlagPSD:end),'g.',[-1.5 1.5],[-1.5 1.5],'k--','LineWidth',2);
    grid;axis([-1.5 1.5 -1.5 1.5]);xlabel('PSD W platmot(m/s)');
    legend('UConn vs PSD','UConnPSD vs PSD','Location','NorthWest')


    %Compares corrected UVW
    XlagU = findlag(UConn_uvw(1,:),PSD_uvwi(1,:),50,0)
    XlagV = findlag(UConn_uvw(2,:),PSD_uvwi(2,:),50,0)
    XlagW = findlag(UConn_uvw(3,:),PSD_uvwi(3,:),50,0)

    figure(4);plot(jdUConn,PSD_uvwi(1,:),'b.--',jdUConn-abs(XlagU/10)/86400,UConn_uvw(1,:),'g*--',jdUConn-abs(XlagU/10)/86400,UConn_uvwPSD(1,:),'m*--',jdUConn-abs(XlagW/10)/86400,UConn_uvwmot(1,:),'k-',jdUConn-abs(XlagW/10)/86400,UConn_uvwmotPSD(1,:),'r-')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('U'' corrected(m/s)');legend('PSD','UConn with GPS','UConnPSD with GPS','UConn without GPS','UConnPSD without GPS')
    figure(5);plot(jdUConn,PSD_uvwi(2,:),'b.--',jdUConn-abs(XlagV/10)/86400,UConn_uvw(2,:),'g*--',jdUConn-abs(XlagV/10)/86400,UConn_uvwPSD(2,:),'m*--',jdUConn-abs(XlagW/10)/86400,UConn_uvwmot(2,:),'k-',jdUConn-abs(XlagW/10)/86400,UConn_uvwmotPSD(2,:),'r-')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('V'' corrected(m/s)');legend('PSD','UConn with GPS','UConnPSD with GPS','UConn without GPS','UConnPSD without GPS')
    figure(6);plot(jdUConn,PSD_uvwi(3,:),'b.--',jdUConn-abs(XlagW/10)/86400,UConn_uvw(3,:),'g*--',jdUConn-abs(XlagW/10)/86400,UConn_uvwPSD(3,:),'m*--',jdUConn-abs(XlagW/10)/86400,UConn_uvwmot(3,:),'k-',jdUConn-abs(XlagW/10)/86400,UConn_uvwmotPSD(3,:),'r-')
    grid;%xlim([jdUConn(1001) jdUConn(2000)])
    ylabel('W'' corrected(m/s)');legend('PSD','UConn with GPS','UConnPSD with GPS','UConn without GPS','UConnPSD without GPS')

    [nanmedian(PSD_uvwi(1,:)) nanmedian(UConn_uvw(1,:)) nanmedian(UConn_uvwPSD(1,:))]
    [nanmedian(PSD_uvwi(2,:)) nanmedian(UConn_uvw(2,:)) nanmedian(UConn_uvwPSD(2,:))]
    [nanmedian(PSD_uvwi(3,:)) nanmedian(UConn_uvw(3,:)) nanmedian(UConn_uvwPSD(3,:))]
end
%Differences in the mean Vj are due to the different methods used. Jim
%compute his fluxes at 20mins and usd high resolution GPS on that window. 
%I use 10min data to correct uvw,and the 10-min means are 0 for Vj
%and Wj. However if I do the average at 20mins, things change a bit (-0.15
%for 328/12 for example).

%computes PSD and CSD for PSD and UConn
[Su,Fu]=psd2(detrend(PSD_uvw(~isnan(PSD_uvw(:,1)),1)),length(PSD_uvw(~isnan(PSD_uvw(:,1)),1)),10,hamming(length(PSD_uvw(~isnan(PSD_uvw(:,1)),1))));[Su,Fu]=specsmoo(Su,10);
[Sv,Fv]=psd2(detrend(PSD_uvw(~isnan(PSD_uvw(:,2)),2)),length(PSD_uvw(~isnan(PSD_uvw(:,2)),2)),10,hamming(length(PSD_uvw(~isnan(PSD_uvw(:,2)),2))));[Sv,Fv]=specsmoo(Sv,10);
[Sw,Fw]=psd2(detrend(PSD_uvw(~isnan(PSD_uvw(:,3)),3)),length(PSD_uvw(~isnan(PSD_uvw(:,3)),3)),10,hamming(length(PSD_uvw(~isnan(PSD_uvw(:,3)),3))));[Sw,Fw]=specsmoo(Sw,10);
[Cuw,Fson_uw]=csd2(detrend(PSD_uvw(~isnan(PSD_uvw(:,3)),3)),detrend(PSD_uvw(~isnan(PSD_uvw(:,1)),1)),length(PSD_uvw(~isnan(PSD_uvw(:,1)),1)),10,hamming(length(PSD_uvw(~isnan(PSD_uvw(:,1)),1))));[Cuw,Fson_uw]=specsmoo(Cuw,10);
[Cvw,Fson_vw]=csd2(detrend(PSD_uvw(~isnan(PSD_uvw(:,3)),3)),detrend(PSD_uvw(~isnan(PSD_uvw(:,2)),2)),length(PSD_uvw(~isnan(PSD_uvw(:,1)),1)),10,hamming(length(PSD_uvw(~isnan(PSD_uvw(:,1)),1))));[Cvw,Fson_vw]=specsmoo(Cvw,10);

[SuUC,FuUC]=psd2(detrend(UConn_uvw(1,:)),length(UConn_uvw(1,:)),10,hamming(length(UConn_uvw(1,:))));[SuUC,FuUC]=specsmoo(SuUC,10);
[SvUC,FvUC]=psd2(detrend(UConn_uvw(2,:)),length(UConn_uvw(2,:)),10,hamming(length(UConn_uvw(2,:))));[SvUC,FvUC]=specsmoo(SvUC,10);
[SwUC,FwUC]=psd2(detrend(UConn_uvw(3,:)),length(UConn_uvw(3,:)),10,hamming(length(UConn_uvw(3,:))));[SwUC,FwUC]=specsmoo(SwUC,10);
[CuwUC,Fson_uwUC]=csd2(detrend(UConn_uvw(3,:)),detrend(UConn_uvw(1,:)),length(UConn_uvw(1,:)),10,hamming(length(UConn_uvw(1,:))));[CuwUC,Fson_uwUC]=specsmoo(CuwUC,10);
[CvwUC,Fson_vwUC]=csd2(detrend(UConn_uvw(3,:)),detrend(UConn_uvw(2,:)),length(UConn_uvw(1,:)),10,hamming(length(UConn_uvw(1,:))));[CvwUC,Fson_vwUC]=specsmoo(CvwUC,10);

% specview(ii,:)=[Fu Su Sv Sw Cuw Cvw FuUC SuUC SvUC SwUC CuwUC CvwUC];
if length(Fu)<54
    Su(54)=NaN;Sv(54)=NaN;Sw(54)=NaN;
    Fu(54)=NaN;Fv(54)=NaN;Fw(54)=NaN;Cuw(54)=NaN;Cvw(54)=NaN;
end
Suview(ii,:)=[Fu Su Sv Sw Cuw Cvw FuUC SuUC SvUC SwUC CuwUC CvwUC];

% figure;clf;semilogx(Fson_uw,Fson_uw.*Cuw,'b',Fson_vw,Fson_vw.*Cvw,'b--',Fson_uwUC,Fson_uwUC.*CuwUC,'r',Fson_vwUC,Fson_vwUC.*CvwUC,'r--',[8e-4 6],[0 0]);grid;xlim([8e-4 6])
figure(1);loglog(Fu,Fu.*Su,'b',FuUC,FuUC.*SuUC,'r')
grid;ylabel('Su(f)');title('U'' PSD');legend('PSD','UConn')
figure(2);clf;loglog(Fv,Fv.*Sv,'b',FvUC,FvUC.*SvUC,'r')
grid;ylabel('Sv(f)');title('V'' PSD');legend('PSD','UConn')
figure(3);clf;loglog(Fw,Fw.*Sw,'b',FwUC,FwUC.*SwUC,'r')
grid;ylabel('Sw(f)');title('W'' PSD');legend('PSD','UConn')

end


Fu=nanmedian(Suview(:,1:54))
Su 
Sv 
Sw 
Cuw 
Cvw 
FuUC 
SuUC 
SvUC 
SwUC 
CuwUC 
CvwUC

lll=54;
DOY=SA(:,1);
fson=SA(:,2:lll+1);
Su=SA(:,lll+2:lll+1+48);
Sv=SA(:,lll*2+2:lll*2+1+48);
Sw=SA(:,lll*3+2:lll*3+1+48);
St=SA(:,lll*4+2:lll*4+1+48);
So=SA(:,lll*5+2:lll*5+1+48);
Suj=SA(:,lll*6+2:lll*6+1+48);
Svj=SA(:,lll*7+2:lll*7+1+48);
Swj=SA(:,lll*8+2:lll*8+1+48);
Cuw=SA(:,lll*9+2:lll*9+1+48);
Cvw=SA(:,lll*10+2:lll*10+1+48);
Ctw=SA(:,lll*11+2:lll*11+1+48);
Cow=SA(:,lll*12+2:lll*12+1+48);
Cuwj=SA(:,lll*13+2:lll*13+1+48);
Cvwj=SA(:,lll*14+2:lll*14+1+48);
