close('all');
fclose('all');
clear;

datasetv='_v1';
read_parameters_HiWinGS_2013;

expp='HiWinGS_2013';
drive='C:\data\CRUISES\HiWinGS_2013\Knorr\flux\Raw';
% drive_w='C:\Data\CRUISES\DYNAMO_2011\flux_data';  
drive_w='C:\data\CRUISES\HiWinGS_2013\flux_data';  

graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device


clearit=1;
printit=1;
year='2013';

yyyy='2013';

decorr=1;  %to decorrelate winds from motion...
dcname=[];
if decorr==1
    dcname='_vdecorr';
end;

doy_PSD=[];accx_PSD=[];accy_PSD=[];accz_PSD=[];ratex_PSD=[];ratey_PSD=[];ratez_PSD=[];
doy_UH=[];accx_UH=[];accy_UH=[];accz_UH=[];ratex_UH=[];ratey_UH=[];ratez_UH=[];
Xlag_all=[];Xlag_all_PML=[];
doy_GPS=[]; pitch_GPS=[];
doy_SCS=[];Hor_acc_SCS=[];Vert_acc_SCS=[];Heave_SCS=[];roll_SCS=[];pitch_SCS=[];
paccx_all = [];paccy_all = [];paccz_all = [];
pratex_all = [];pratey_all = [];pratez_all = [];
doy_PML=[];accx_PML=[];accy_PML=[];accz_PML=[];ratex_PML=[];ratey_PML=[];ratez_PML=[];
All_accx_PSDPML=[];All_accx_PML=[];All_accx_UH=[];All_accx_PSDUH=[];
All_accy_PSDPML=[];All_accy_PML=[];All_accy_UH=[];All_accy_PSDUH=[];
All_accz_PSDPML=[];All_accz_PML=[];All_accz_UH=[];All_accz_PSDUH=[];

rge=283;
for ddd=rge
    [m,d]=yd2md(str2num(year), ddd);
    Vdate=[str2num(year) m d];
%     jd=sprintf('%03i',ddd);
   
    for hhh=00;   %
%        hr=sprintf('%02i',hhh);
       read_motion_HiWinGS_2013;  %
       [jd_UH,AccX_UH,AccY_UH,AccZ_UH,RateX_UH,RateY_UH,RateZ_UH,U_UH,V_UH,W_UH,Tsonic_UH,LagPulse,Xlag]=read_dms_HiWinGS_2013(ddd,hhh,yyyy,accz);

       doy_PSD=[doy_PSD ddd+(hhh+(tref+Xlag/10)/3600)/24];
       accx_PSD=[accx_PSD accx];accy_PSD=[accy_PSD accy];accz_PSD=[accz_PSD accz];
       ratex_PSD=[ratex_PSD ratex];ratey_PSD=[ratey_PSD ratey];ratez_PSD=[ratez_PSD ratez];

       doy_UH=[doy_UH jd_UH];
       accx_UH=[accx_UH AccX_UH'];accy_UH=[accy_UH AccY_UH'];accz_UH=[accz_UH AccZ_UH'];
       ratex_UH=[ratex_UH RateX_UH'];ratey_UH=[ratey_UH RateY_UH'];ratez_UH=[ratez_UH RateZ_UH'];

       Xlag_all=[Xlag_all Xlag]; 
        All_accx_UH=[All_accx_UH; AccX_UH(1:length(accx)-abs(Xlag))];
        All_accx_PSDUH=[All_accx_PSDUH accx(1+abs(Xlag):length(accx))];     
        All_accy_UH=[All_accy_UH; AccY_UH(1:length(accy)-abs(Xlag))];
        All_accy_PSDUH=[All_accy_PSDUH accy(1+abs(Xlag):length(accy))];     
        All_accz_UH=[All_accz_UH; AccZ_UH(1:length(accz)-abs(Xlag))];
        All_accz_PSDUH=[All_accz_PSDUH accz(1+abs(Xlag):length(accz))];     
        
       read_head_pitch_HiWinGS_2013; % 
       doy_GPS=[doy_GPS jd_hed_pitch'];
       pitch_GPS=[pitch_GPS pitchx_pitch];

       [jd_gps,Hor_acc,Vert_acc,Heave,Roll,Pitch]=read_rawSCS_MRU(ddd,hhh,year);     
       doy_SCS=[doy_SCS jd_gps];
       Hor_acc_SCS=[Hor_acc_SCS Hor_acc];
       Vert_acc_SCS=[Vert_acc_SCS Vert_acc];
       Heave_SCS=[Heave_SCS Heave];
       roll_SCS=[roll_SCS Roll];
       pitch_SCS=[pitch_SCS -Pitch];

        paccx_all = [paccx_all; polyfit(accx(1+abs(Xlag):end),AccX_UH(1:end-abs(Xlag))',1)];
        paccy_all = [paccy_all; polyfit(accy(1+abs(Xlag):end),AccY_UH(1:end-abs(Xlag))',1)];
        paccz_all = [paccz_all; polyfit(accz(1+abs(Xlag):end),AccZ_UH(1:end-abs(Xlag))',1)];
        pratex_all = [pratex_all; polyfit(ratex(1+abs(Xlag):end),RateX_UH(1:end-abs(Xlag))',1)];
        pratey_all = [pratey_all; polyfit(ratey(1+abs(Xlag):end),RateY_UH(1:end-abs(Xlag))',1)];
        pratez_all = [pratez_all; polyfit(ratez(1+abs(Xlag):end),RateZ_UH(1:end-abs(Xlag))',1)];
        
        %adds in PML data
        [jd_PML,AccX_PML,AccY_PML,AccZ_PML,RateX_PML,RateY_PML,RateZ_PML,U_PML,V_PML,W_PML,Tsonic_PML,Xlag]=read_PML_HiWinGS_2013(ddd,hhh,yyyy,accz);
        doy_PML=[doy_PML jd_PML-Xlag/10/3600/24];
        accx_PML=[accx_PML AccX_PML];accy_PML=[accy_PML AccY_PML];accz_PML=[accz_PML AccZ_PML];
        ratex_PML=[ratex_PML RateX_PML];ratey_PML=[ratey_PML RateY_PML];ratez_PML=[ratez_PML RateZ_PML];

        Xlag_all_PML=[Xlag_all_PML Xlag];
        
        if length(AccZ_PML)>36000
            All_accx_PSDPML=[All_accx_PSDPML  accx(1+abs(Xlag):length(accx))];
            All_accx_PML=[All_accx_PML  AccX_PML(1:length(accx)-abs(Xlag))];
            All_accy_PSDPML=[All_accy_PSDPML  accy(1+abs(Xlag):length(accy))];
            All_accy_PML=[All_accy_PML  AccY_PML(1:length(accy)-abs(Xlag))];
            All_accz_PSDPML=[All_accz_PSDPML  accz(1+abs(Xlag):length(accz))];
            All_accz_PML=[All_accz_PML  AccZ_PML(1:length(accz)-abs(Xlag))];
        else
            All_accx_PSDPML=[All_accx_PSDPML  accx(1+abs(Xlag):(length(AccX_PML)))];
            All_accx_PML=[All_accx_PML  AccX_PML(1:end-abs(Xlag))];      
            All_accy_PSDPML=[All_accy_PSDPML  accy(1+abs(Xlag):(length(AccY_PML)))];
            All_accy_PML=[All_accy_PML  AccY_PML(1:end-abs(Xlag))];            
            All_accz_PSDPML=[All_accz_PSDPML  accz(1+abs(Xlag):(length(AccZ_PML)))];
            All_accz_PML=[All_accz_PML  AccZ_PML(1:end-abs(Xlag))];            
        end    
    end
end

doy_UH(doy_UH<282.5)=NaN;
accz_PSD(accz_PSD<0)=NaN;
% ratex_PSD(ratex_PSD>0.2)=NaN;ratex_PSD(ratex_PSD<-1)=NaN;
ratey_PSD(ratey_PSD>0.2)=NaN;ratey_PSD(ratey_PSD<-1)=NaN;
ratez_PSD(ratez_PSD>0.2)=NaN;

figure;plot(doy_PSD,accz_PSD,doy_UH,accz_UH,'r',doy_PML,accz_PML,'k');title('accz');
figure;plot(accz_PSD(1+abs(Xlag_all(1)):end),accz_UH(1:end-abs(Xlag_all(1))),'.')
% figure;plot(accz_PSD(1+abs(Xlag_all(2))+36000:end),accz_UH(1+36000:end-abs(Xlag_all(2))),'.')
% figure;plot(accz_PSD(36000:end),accz_UH(36000:end),'.');title('accz');
% axis([-0.1 0.1 -0.1 0.1])
%Note: Interestingly the first hour is the one with most scatter in the
%VS plot...need to check why 
% figure;plot(ratey_PSD(1:36000),ratey_UH(1:36000),'.');title('accz');
% --> it's because of the lag not taken into account!
% figure;plot(accz_PSD(1+abs(Xlag_all(1)):36000),accz_UH(1:36000-abs(Xlag_all(1))),'.')

% 
% for hhh=00:01
% %     figure(hhh+1);plot(accz_PSD(36000*hhh+1:36000*(hhh+1)),accz_UH(36000*hhh+1:36000*(hhh+1)),'.');title('accz');
%     figure;plot(accz_PSD(36000*hhh+1+abs(Xlag_all(hhh+1)):36000*(hhh+1)),accz_UH(36000*hhh+1:36000*(hhh+1)-abs(Xlag_all(hhh+1))),'.')
% hold on;plot([0 16],[0 16],'k--')
% end


% hhh=0;figure;plot(doy_PSD(36000*hhh+1:36000*(hhh+1)),accz_PSD(36000*hhh+1:36000*(hhh+1)),doy_UH(36000*hhh+1:36000*(hhh+1)),accz_UH(36000*hhh+1:36000*(hhh+1)),'r');title('accz'); 

[nanmedian(accx_UH) nanmedian(accx_PSD) nanmedian(accx_PML) nanmedian(accy_UH) nanmedian(accy_PSD) nanmedian(accy_PML) nanmedian(accz_UH) nanmedian(accz_PSD) nanmedian(accz_PML)]
[nanmedian(ratex_UH) nanmedian(ratex_PSD) nanmedian(ratex_PML) nanmedian(ratey_UH) nanmedian(ratey_PSD) nanmedian(ratey_PML) nanmedian(ratez_UH) nanmedian(ratez_PSD) nanmedian(ratez_PML)]
[nanmedian(accx_UH-accx_PSD) nanmedian(accy_UH-accy_PSD) nanmedian(accz_UH-accz_PSD)]
[nanmedian(ratex_UH-ratex_PSD) nanmedian(ratey_UH-ratey_PSD) nanmedian(ratez_UH-ratez_PSD)]

G_PSD=sqrt(nanmean(accx_PSD)^2+nanmean(accy_PSD)^2+nanmean(accz_PSD)^2);
G_UH=sqrt(nanmean(accx_UH)^2+nanmean(accy_UH)^2+nanmean(accz_UH)^2);
G_PML=sqrt(nanmean(accx_PML)^2+nanmean(accy_PML)^2+nanmean(accz_PML)^2);


A =importdata('S:\cruises\HIWINGS_2013\Knorr\Scientific_analysis\MeanNBulkfluxes\HiWinGS_2013_hrMeansNFluxes_r1.txt');
A=A.data;
yday=A(:,1);    %Decimal yearday (UTC)
Lat=A(:,2);     %Latitude (deg)
SOG=A(:,4);     %Speed over ground (m/s)
G_lat=nanmean(grv(Lat(yday>=doy_PSD(1) & yday<doy_PSD(end))))

[G_PSD G_UH G_lat G_PML]

% figure;plot(yday,SOG,'+--');xlim([ddd ddd+1])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Look at the tilts while in port
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bf=ones(1,80)/80;          %incr/2 point averaging filter
af=1;

roll_PSD=atan2(filtfilt(bf,af,accy_PSD),G_lat)*180/pi;
pitch_PSD = atan2(filtfilt(bf,af,-accx_PSD),G_lat)*180/pi;
% roll_PSD2=cumtrapz(detrend(ratex_PSD))*0.05;
% pitch_PSD2=cumtrapz(detrend(ratey_PSD))*0.05;

roll_UH=accy_UH;pitch_UH=accy_UH;
roll_UH(~isnan(accy_UH))=atan2(filtfilt(bf,af,accy_UH(~isnan(accy_UH))),G_lat)*180/pi;
pitch_UH(~isnan(accx_UH)) = atan2(filtfilt(bf,af,-accx_UH(~isnan(accx_UH))),G_lat)*180/pi;
% roll_UH2=accy_UH;pitch_UH2=accy_UH;
% roll_UH2(~isnan(ratex_UH))=cumtrapz(detrend(ratex_UH(~isnan(ratex_UH))))*0.05;
% pitch_UH2(~isnan(ratey_UH))=cumtrapz(detrend(ratey_UH(~isnan(ratey_UH))))*0.05;
% 
roll_PML=accy_PML;pitch_PML=accy_PML;
roll_PML(~isnan(accy_PML))=atan2(filtfilt(bf,af,accy_PML(~isnan(accy_PML))),G_lat)*180/pi;
pitch_PML(~isnan(accx_PML)) = atan2(filtfilt(bf,af,-accx_PML(~isnan(accx_PML))),G_lat)*180/pi;

figure;plot(doy_SCS,filtfilt(bf,af,roll_SCS),'k',doy_PSD,roll_PSD,doy_UH,roll_UH,'r',doy_PML,roll_PML,'g');title('roll (deg)'); legend('SCS','PSD','UH','PML');ylim([-7 7])
roll_SCSf=filtfilt(bf,af,roll_SCS);
% [nanmedian(roll_PSD(doy_PSD<282.8)) nanmedian(roll_UH(doy_UH<282.8)) nanmedian(roll_SCSf(doy_SCS<282.8))]
[nanmedian(roll_PSD) nanmedian(roll_UH) nanmedian(roll_PML) nanmedian(roll_SCSf) ]

pitch_SCSf=filtfilt(bf,af,pitch_SCS);pitch_SCSf(pitch_SCSf>0.5)=NaN;
pitch_GPSf=pitch_GPS;pitch_GPSf(~isnan(pitch_GPS))=filtfilt(bf,af,pitch_GPS(~isnan(pitch_GPS)));
figure;plot(doy_SCS,filtfilt(bf,af,pitch_SCS),'k',doy_PSD,pitch_PSD,doy_UH,pitch_UH,'r',doy_GPS(~isnan(pitch_GPSf)),pitch_GPSf(~isnan(pitch_GPSf)),'m',doy_PML,pitch_PML,'g');title('pitch (deg)'); legend('SCS','PSD','UH','diffGPS','PML')
% [nanmedian(pitch_PSD(doy_PSD<282.8)) nanmedian(pitch_UH(doy_UH<282.8)) nanmedian(pitch_GPS(doy_GPS<282.8)) nanmedian(pitch_SCSf(doy_SCS<282.8))]
[nanmedian(pitch_PSD) nanmedian(pitch_UH) nanmedian(pitch_PML) nanmedian(pitch_GPS) nanmedian(pitch_SCSf)]


figure;plot(accz_PSD(1+abs(Xlag_all(1)):end),accz_UH(1:end-abs(Xlag_all(1))),'.')

figure;
subplot(3,2,1);plot(doy_PSD,detrend(accx_PSD),doy_UH,detrend(accx_UH),'r',doy_PML+4/10/3600/24,detrend(accx_PML),'k');title('accx');xlim([284.3697 284.3706]);grid
subplot(3,2,2);plot(All_accx_PSDPML,All_accx_PML,'c.',All_accx_PSDUH,All_accx_UH,'.');hold on;plot([-5 5],[-5 5],'k--');axis([-5 5 -5 5]);grid
subplot(3,2,3);plot(doy_PSD,detrend(accy_PSD),doy_UH,detrend(accy_UH),'r',doy_PML+4/10/3600/24,detrend(accy_PML),'k');title('accy');xlim([284.3697 284.3706]);grid
subplot(3,2,4);plot(All_accy_PSDPML,All_accy_PML,'c.',All_accy_PSDUH,All_accy_UH,'.');hold on;plot([-5 5],[-5 5],'k--');axis([-5 5 -5 5]);grid
% subplot(3,3,7);plot(doy_PSD,detrend(accz_PSD),doy_UH,detrend(accz_UH),'r',doy_PML+4/10/3600/24,detrend(accz_PML),'k');title('accz');xlim([ddd ddd+1]);grid
subplot(3,2,5);plot(doy_PSD,detrend(accz_PSD),doy_UH,detrend(accz_UH),'r',doy_PML+4/10/3600/24,detrend(accz_PML),'k');title('accz');xlim([284.3697 284.3706]);grid
subplot(3,2,6);plot(All_accz_PSDPML,All_accz_PML,'c.',All_accz_PSDUH,All_accz_UH,'.');hold on;plot([0 20],[0 20],'k--');axis([0 20 0 20]);grid
