% compare ship and PSD measurements for the cruise

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSDmet=load('G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_met_1min_all.txt');
jd_pc=PSDmet(:,1);  %decimal JD
PSP=PSDmet(:,3);
PIR=PSDmet(:,2);
Tc=PSDmet(:,4);
Td=PSDmet(:,5);
Tsnk=PSDmet(:,6);
Tvais=PSDmet(:,7);    %degC
Rhvais=PSDmet(:,8);    % %
press=PSDmet(:,12);    % mbar
org=PSDmet(:,9);    % mbar
Qvais=qair3([press Tvais Rhvais]); % min(  ,100) rejects supersaturation
Qvais(Qvais>50)=NaN;
radmet=load('G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_rad_1min_all.txt');
Tc1=radmet(:,6);
Td1=radmet(:,7);
Tc2=radmet(:,8);
Td2=radmet(:,9);


SCSmet=load('G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_scs_1min_all.txt');
[c, ia, ib] = intersect(jd_pc,SCSmet(:,1));  %indexes common to both AB and AC
% SCSmet=SCSmet(1:end-1,:);
WXTP_PSP=SCSmet(:,7);WXTS_PIR=SCSmet(:,8);
WXTP_Ta=SCSmet(:,9);WXTS_Ta=SCSmet(:,10);
WXTP_Pa=SCSmet(:,12);WXTS_Pa=SCSmet(:,11);
WXTP_Ua=SCSmet(:,22);WXTS_Ua=SCSmet(:,21);
QscsP=qair3([WXTP_Pa WXTP_Ta WXTP_Ua]); % min(  ,100) rejects supersaturation
QscsS=qair3([WXTS_Pa WXTS_Ta WXTS_Ua]); % min(  ,100) rejects supersaturation
WXTP_TS=SCSmet(:,27);WXTS_TS=SCSmet(:,28);
WXTP_TD=SCSmet(:,25);WXTS_TD=SCSmet(:,26);
WXTP_sst=SCSmet(:,6);
WXTP_lat=SCSmet(:,2);WXTP_lon=SCSmet(:,3);
WXTP_cog=SCSmet(:,4);WXTP_sog=SCSmet(:,5)*0.514;
WXTP_head=SCSmet(:,23);
WXTP_Sm=SCSmet(:,20);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Look at pressure
%PSD height is 11m from WL
%SCS assumes 15m from WL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zSCS=15;zPSD=11;
WXTP_Pa(WXTP_Pa<800)=NaN;WXTS_Pa(WXTS_Pa<800)=NaN;press(press<800)=NaN;
figure;plot(jd_pc,press+0.12*(zPSD),'g',jd_pc,WXTP_Pa+0.12*(zSCS),'b',jd_pc,WXTS_Pa+0.12*(zSCS),'r');grid;
ylim([950 1050])

figure;plot(jd_pc,WXTP_Pa+0.12*(zSCS)-(press+0.12*(zPSD)),'b',jd_pc,WXTS_Pa+0.12*(zSCS)-(press+0.12*(zPSD)),'r');grid;
[nanmedian(WXTP_Pa+0.12*(zSCS)-(press+0.12*(zPSD))) nanmedian(WXTS_Pa+0.12*(zSCS)-(press+0.12*(zPSD)))]

figure;plot(WXTS_TS,WXTP_Pa+0.12*(zSCS)-(press+0.12*(zPSD)),'b.',WXTS_TS,WXTS_Pa+0.12*(zSCS)-(press+0.12*(zPSD)),'r.');grid;
legend('Port','Stbd');xlabel('WS')

figure;subplot(2,1,1);plot(jd_pc,press+0.12*(zPSD),'g',jd_pc,WXTP_Pa+0.12*(zSCS),'b',jd_pc,WXTS_Pa+0.12*(zSCS),'r');grid;
xlim([305 306])
title('sea level BP');legend('PSD','SCS port','SCS stbd')
subplot(2,1,2);plot(jd_pc,WXTS_TS);grid;xlim([305 306]);title('WS')
%PSD seems to be affected by dynamic pressure fluctuations...will use SCS measurement
figure;plot(WXTP_Pa,WXTS_Pa,'.');grid;xlabel('port');ylabel('stbd')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Look at SST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%next, clean up times when snkae off the air and replaces values with SCS
%data -0.1degC
% Tsnk(jd_pc<283.1)=WXTP_sst(jd_pc<283.1)-0.1;
% Tsnk(jd_pc>288.1 & jd_pc<288.7)=WXTP_sst(jd_pc>288.1 & jd_pc<288.7)-0.1;
% Tsnk(jd_pc>298.06 & jd_pc<298.09)=WXTP_sst(jd_pc>298.06 & jd_pc<298.09)-0.1;
% Tsnk(jd_pc>298.885 & jd_pc<298.9)=WXTP_sst(jd_pc>298.885 & jd_pc<298.9)-0.3;
Tsnk(jd_pc<283.4)=NaN;
Tsnk(jd_pc>288.1 & jd_pc<288.7)=NaN;
Tsnk(jd_pc>298.06 & jd_pc<298.09)=NaN;
Tsnk(jd_pc>298.885 & jd_pc<298.9)=NaN;
WXTP_sst(jd_pc>306.6 & jd_pc<306.7)=NaN;
Tsnk(Tsnk>20)=NaN;

Tsnk=Tsnk+0.096;
PSDtsnkf=Tsnk;PSDtsnkf(~isnan(Tsnk))=filtfilt(bf,af,Tsnk(~isnan(Tsnk)));

figure;plot(jd_pc,WXTS_TS,'k--',jd_pc,Tsnk,'g',jd_pc,WXTP_sst,'r',jd_pc,WXTS_Ta);grid;
legend('SCS WS','PSD snake','SCS Tsea','SCS Tair')
xlabel('DOY (2013)')

figure;plot(WXTS_TS,Tsnk-WXTP_sst,'b.');grid;
xlabel('WS')
figure;subplot(2,1,1);plot(jd_pc,Tsnk,'g',jd_pc,WXTP_sst,'r');grid;
% xlim([298 300])
title('Tsea');legend('PSD (0.05m','SCS (5m)')
subplot(2,1,2);plot(jd_pc,WXTS_TS);grid;title('WS')
% xlim([298 300]);
%Snake and SCS within 0.1degC when snake in the water -->will use SST sea
%temperature! but some features seems differente too (cf 291.8) so might
%correct snake for wind issues.
y=(Tsnk-WXTP_sst);x=WXTS_TS;binsize=0.5;xmin=0;xmax=32;
figure;plot(x,y,'b.');grid;
xlabel('WS')
% y(y<-1)=NaN;y(y>.2)=NaN;
% y=y+0.1;
% hold on;plot(x,y,'g.')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(x,y,binsize,xmin,xmax);
hold on;plot(X,Y,'rd--')

% Tsnk-WXTP_sst = 5.8e-6*x.^3 - 0.00088*x.^2 + 0.013*x -0.1
% Tsnk-WXTP_sst = -6e-4*x.^2 + 0.0095*x -0.13
% with 0.1 offset off
% Tsnk-WXTP_sst = 5.8e-6*x.^3 - 0.00088*x.^2 + 0.013*x -0.038
% Tsnk-WXTP_sst = -6e-4*x.^2 + 0.0095*x -0.029
% --> I like hold on;plot(jd_pc,Tsnk-(-6e-4*x.^2 + 0.0095*x -0.029 ),'y')
%after removing snake out of water events and without 0.1C offset
% Tsnk-WXTP_sst = -8.7e-6*x.^3 - 0.00022*x.^2 + 0.0061*x -0.022
% Tsnk-WXTP_sst = -6.3e-4*x.^2 + 0.011*x -0.036
% Tsnk-WXTP_sst = -4.6e-4*x.^2 + 0.0066*x -0.015
%do some smoothing
N=5;
bf=ones(1,N)/N;  %Use this for 1 minute output
af=1;
% PSDtsnkf=filtfilt(bf,af,Tsnk-(-6e-4*x.^2 + 0.0095*x -0.029 ));

% PSDtsnkf=cleanup_lb(Tsnk-(-6e-4*x.^2 + 0.0095*x -0.029 ),jd_pc,jd_pc,find(~isnan(Tsnk)),5);
PSDtsnkf=cleanup_lb(Tsnk-(-6.3e-4*x.^2 + 0.011*x -0.036),jd_pc,jd_pc,find(~isnan(Tsnk)),5);

hold on;plot(jd_pc,PSDtsnkf,'m')

filesave='C:\data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\cleanedTsea';
eval(['save ' filesave ' jd_pc Tsnk WXTP_sst'])  
clear filesave

%do the same against detT. Around doy288 snake higher as Tair warmer -->
%might be slightly affected by flying time...
y=(Tsnk-WXTP_sst);
x=WXTP_sst-WXTS_Ta;
binsize=1;xmin=-4;xmax=15;
figure;plot(x,y,'b.');grid;
xlabel('Tsea-Tair')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(x,y,binsize,xmin,xmax);
hold on;plot(X,Y,'rd--')
% Tsnk-WXTP_sst = -2.4e-4*x.^2 - 0.0028*x + 0.0017
% Tsnk-WXTP_sst = -0.0054*x + 0.0018
% PSDtsnkf2=Tsnk-(-2.4e-4*x.^2 - 0.0028*x + 0.0017);
PSDtsnkf2=Tsnk-(-0.0054*x + 0.0018);


y=(Tsnk-WXTP_sst);x=WXTS_TS;detT=WXTP_sst-WXTS_Ta;
figure;scatter(x,y,3,detT);grid;
h = colorbar;ylabel(h, 'detT')
xlabel('WS');ylabel('Tsnk - SCS Tsea');

y=(Tsnk-WXTP_sst);x=WXTP_sst-WXTS_Ta;
figure;scatter(x,y,3,WXTP_TS);grid;
h = colorbar;ylabel(h, 'WS')
xlabel('SCS Tsea-Tair');ylabel('Tsnk - SCS Tsea');
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(x,y,1,-5,15);
hold on;plot(X,Y,'rd--')


x=WXTP_sst-WXTS_Ta;
y=(Tsnk-WXTP_sst);
%averaged by detT bins
figure
uu=-5:2:15;
for i=1:(length(uu)-1)
  ii=find(x>=uu(i) & x<uu(i+1));
  [X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(WXTS_TS(ii),y(ii),1,0,30);
  col=rand(3,1);
  hold on;plot(X,Y,'Color',col);hold on;
  legendInfo{i} = sprintf('detT = %i - %i',uu(i),uu(i+1));
end; 
[legend_h,object_h,plot_h,text_strings] = legend(legendInfo);


% Hs_heave2wavespline_all=load('G:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_waveplineinterpolatedtoheave_minusriegldplusheave_vcfrp\Hs_hourly_all.prn');
% doy_all=Hs_heave2wavespline_all(:,1)+Hs_heave2wavespline_all(:,2)/24;
% ii=find(abs(Hs_heave2wavespline_all(:,4))<7 );
% figure;plot(doy_all(ii),Hs_heave2wavespline_all(ii,4))


figure;plot(PSP,Tsnk-WXTP_sst,'r.')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compares lat, lon, sog, cog, heading
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSDnav=load('G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_gpsnav_1min_all.txt');
jd_nav=PSDnav(:,1);
PSDlat=PSDnav(:,2);%decimal latitude, deg
PSDlon=PSDnav(:,3);%decimal longitude, deg
PSDsog=PSDnav(:,4);%GPS SOG, m/s
PSDcog=PSDnav(:,5);%GPS COG, deg
PSDhead=PSDnav(:,6);%Crescent GPS heading, deg

PSDhead2=PSDhead;  
PSDhead2(jd_nav<296.7)=PSDhead(jd_nav<296.7)-1.0685;
PSDhead2(jd_nav>300)=PSDhead(jd_nav>300)-0.715;
PSDhead2(jd_nav>298.646 & jd_nav<300.6)=WXTP_head(jd_pc>298.646 & jd_pc<=300.6);


figure;plot(jd_nav,PSDlat,jd_pc,WXTP_lat,'g'); %OK
figure;plot(jd_nav,PSDlon,jd_pc,WXTP_lon,'g'); %OK
figure;plot(jd_pc,WXTP_sog,'g',jd_nav,PSDsog); %SCS sog pretty noisy... will use PSD
figure;plot(jd_pc,WXTP_cog,'g',jd_nav,PSDcog); %SCS cog pretty noisy... will use PSD
figure;plot(jd_pc,WXTP_head,'g',jd_nav,PSDhead2); %PSD and SCS within 1-2deg for heading. Willuse SCS as more aligned hopefully...
figure;plot(jd_pc(1:end-1),WXTP_head(1:end-1)-PSDhead2); %


SPD=SCSmet(:,end)*0.514; %SOW?
CurrSpd=PSDsog-SPD(1:end-1);
CurrSpd2=WXTP_sog-SPD;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compares WS and direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSDson=load('G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_son_1min_all.txt');
jd_son=PSDson(:,1);
PSDu=PSDson(:,2);
PSDv=PSDson(:,3);
PSDw=PSDson(:,4);
PSDdir=PSDson(:,6);

% S=sqrt(u.*u+v.*v+w.*w);
% dir=atan2(-V,U)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
% dir=mod(dir+720,360);
% uu=jazx1(2,:)/.95;				%.95 is empirical distortion factor
% vv=jazx1(3,:)/1.15;				%1.15 is other empirical distortion factor
rdcon=pi/180;
lrgc=cos(PSDhead2*rdcon);
lrgs=sin(PSDhead2*rdcon);
un=PSDu.*lrgc+PSDv.*lrgs;
ue=PSDu.*lrgs-PSDv.*lrgc;
sogx=PSDsog.*cos(PSDcog*rdcon);			%sogx
sogy=PSDsog.*sin(PSDcog*rdcon);			%sogy

step=1/1440;      %1 minute time step
yday=282:step:318;
ii=find(~isnan(sogx));sogx=interp1(jd_nav(ii),sogx(ii),jd_nav);
ii=find(~isnan(sogy));sogy=interp1(jd_nav(ii),sogy(ii),jd_nav);
ii=find(~isnan(lrgc));lrgc=interp1(jd_nav(ii),lrgc(ii),jd_nav);
ii=find(~isnan(lrgs));lrgs=interp1(jd_nav(ii),lrgs(ii),jd_nav);

unt=un+sogx;     
uet=ue+sogy;     
PSDst=sqrt(unt.*unt+uet.*uet);
PSDtruedir=mod(atan2(-uet,-unt)/rdcon+360,360);

% figure;plot(jd_son,PSDst,'g',jd_pc,WXTP_TS,jd_pc,WXTS_TS,'r');grid;
% legend('PSD','SCS port','SCS stbd')
% figure;
% subplot(2,2,1);plot(PSDst,WXTP_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;
% xlabel('PSD 15.9m');ylabel('SCS port 15m')
% subplot(2,2,2);plot(PSDst,WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;
% xlabel('PSD 15.9m');ylabel('SCS stbd 15m')
% subplot(2,2,3);plot(WXTP_TS(1:end-1),WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;
% xlabel('SCS port 15m');ylabel('SCS stbd 15m')

%get BB true wind
% % rdcon=pi/180;
% % for doy=298:317
% %     [y, m, d, ~, ~, ~] = datevec(doy+datenum(2013,1,0));  
% % %     read_UHwind_HiWinGS_2013  %motion corrected data, rotated into stream
% %     read_dms_HiWinGS_2013_daily
% %     
% %     PSDmet=load(['G:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_proc_met_1min_' sprintf('%04i_%02i_%02i_%03i',y,m,d,doy) '.txt']);
% %     ii=find(~isnan(jd_UH));
% %     jd_pc=PSDmet(:,1);  %decimal JD
% %     Ui=interp1(jd_UH(ii),U_UH(ii),jd_pc);
% %     Vi=interp1(jd_UH(ii),V_UH(ii),jd_pc);
% %     Wi=interp1(jd_UH(ii),W_UH(ii),jd_pc);
% %     Tsonici=interp1(jd_UH(ii),Tsonic_UH(ii),jd_pc);
% %     
% % %     Dirx=cos(RelWdir_UH'*rdcon);    %convert to N
% % %     Diry=sin(RelWdir_UH'*rdcon);    %convert to E
% % % %     ii=find(~isnan(jd_UH));
% % %     Dirxi=(intavg(jd_UH,Dirx,jd_pc(ii))); 
% % %     Diryi=(intavg(jd_UH,Diry,jd_pc(ii))); 
% % %     RelWdiri=(mod(atan2(Diryi,Dirxi)/rdcon+360,360));
% % 
% %     diri=atan2(-Vi,Ui)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
% %     diri=mod(diri+720,360);
% % 
% %     prt_jas_UHson_HiWinGS_2013
% %  
% % end

UHson=load('C:\data\CRUISES\HiWinGS_2013\Knorr\UH_DMS\HiWinGS_2013_proc_UHson_1min_all.txt');
jd_UH=UHson(1:end-1,1);
UHu=UHson(1:end-1,2);
UHv=UHson(1:end-1,3);

% UHu=UHson(1:end-1,2)/.95;
% UHv=UHson(1:end-1,3)/1.15;
UHw=UHson(1:end-1,4);
UHdir=UHson(1:end-1,5);
% figure;plot(jd_son,PSDu,jd_UH,UHu,'r')
% [c, ia, ib] = intersect(jd_son,jd_UH);  %indexes common to both AB and AC
UH_un=UHu.*lrgc+UHv.*lrgs;
UH_ue=UHu.*lrgs-UHv.*lrgc;
UH_unt=UH_un+sogx;     
UH_uet=UH_ue+sogy;     
UH_st=sqrt(UH_unt.*UH_unt+UH_uet.*UH_uet);
UH_st(isnan(UH_st))=PSDst(isnan(UH_st));
UHtruedir=mod(atan2(-UH_uet,-UH_unt)/rdcon+360,360);
UHtruedir(isnan(UHtruedir))=PSDtruedir(isnan(UHtruedir));

UH_st(jd_UH>300.8095 & jd_UH<300.81)=NaN;
UH_st(jd_UH>301.5095 & jd_UH<301.51)=NaN;
UH_st(jd_UH>304.8002 & jd_UH<304.801)=NaN;
UH_st(jd_UH>304.857 & jd_UH<304.86)=NaN;
UH_st(jd_UH>313.5406 & jd_UH<313.5414)=NaN;

figure;plot(jd_son,PSDst,'g',jd_pc,WXTP_TS,jd_pc,WXTS_TS,'r',jd_UH,UH_st,'k');grid;
legend('PSD','SCS port','SCS stbd','UH')

figure;
subplot(2,2,1);plot(PSDst,WXTP_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15.9m');ylabel('SCS port 15m')
subplot(2,2,2);plot(PSDst,WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15.9m');ylabel('SCS stbd 15m')
subplot(2,2,3);plot(WXTP_TS(1:end-1),WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('SCS port 15m');ylabel('SCS stbd 15m')
subplot(2,2,4);plot(PSDst,UH_st,'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15.9m');ylabel('UH 16.29m')

figure;plot(jd_son,PSDst-UH_st);grid;
WSdiff=PSDst-UH_st;
ll=find(UHdir>180);
figure;plot(UHdir,WSdiff,'b.',UHdir(ll)-360,WSdiff(ll),'b.',[-90 -90],[-20 20],'r',[90 90],[-20 20],'r');axis([-180 180 -10 10]);grid
xlabel('reldir');ylabel('PSD - UH true speed diff (m/s)')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir,WSdiff,5,0,180);
hold on;plot(X,Y,'rd--')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir(ll)-360,WSdiff(ll),5,-180,0);
hold on;plot(X,Y,'rd--')

figure;subplot(2,1,1);plot(jd_nav,PSDsog);grid
subplot(2,1,2);plot(jd_nav,PSDcog,jd_nav,PSDhead2,'r',jd_UH,UHdir,'k.');legend('COG','HEADING','Relative Wndir');grid
orient tall

%do height adjustements to 15m
zi=600*ones(size(PSDst));
zref_u=15;
zref_t=15;
zref_q=15;
RHscst=relhum([WXTS_Ta QscsS/0.93-0.34 WXTP_Pa]);
PSDc1=coare30vn_ref([PSDst; NaN],15.9,WXTP_Ta,15,RHscst,15,WXTP_Pa+0.12*(15),WXTP_sst,PSP,PIR,WXTP_lat,zi,zref_u,zref_t,zref_q,1);
UHc1=coare30vn_ref([UH_st; NaN],16.29,WXTP_Ta,15,RHscst,15,WXTP_Pa+0.12*(15),WXTP_sst,PSP,PIR,WXTP_lat,zi,zref_u,zref_t,zref_q,1);
WSdiff=PSDc1(1:end-1,19)-UHc1(1:end-1,19);
ll=find(UHdir>180);
figure;plot(UHdir,WSdiff,'b.',UHdir(ll)-360,WSdiff(ll),'b.',[-90 -90],[-20 20],'r',[90 90],[-20 20],'r');axis([-180 180 -10 10]);grid
xlabel('reldir');ylabel('PSD - UH true speed diff (m/s)')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir,WSdiff,5,0,180);
hold on;plot(X,Y,'rd--')
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir(ll)-360,WSdiff(ll),5,-180,0);
hold on;plot(X,Y,'rd--')

figure;
subplot(2,2,1);plot(PSDc1(1:end-1,19),WXTP_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15.9m');ylabel('SCS port 15m')
subplot(2,2,2);plot(PSDc1(1:end-1,19),WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15.9m');ylabel('SCS stbd 15m')
subplot(2,2,3);plot(WXTP_TS(1:end-1),WXTS_TS(1:end-1),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('SCS port 15m');ylabel('SCS stbd 15m')
subplot(2,2,4);plot(PSDc1(1:end-1,19),UHc1(1:end-1,19),'.',[0 35],[0 35],'k--');grid;axis([0 35 0 35])
xlabel('PSD 15m');ylabel('UH 15m')


WSdiff=PSDc1(1:end-1,19)-UHc1(1:end-1,19);
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir,WSdiff,5,0,180);
figure;
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
plot(X,Y,'gd--');axis([-180 180 -10 10]);grid
WSdiff=WXTS_TS-UHc1(:,19);
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb([UHdir; NaN],WSdiff,5,0,180);hold on;plot(X,Y,'rd--')
WSdiff=WXTP_TS-UHc1(:,19);
[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb([UHdir; NaN],WSdiff,5,0,180);hold on;plot(X,Y,'bd--')
legend('PSD - UH','SCS stbd - UH','SCS port - UH')
hold on;plot([-180 180],[0 0],'k--')
hold on;plot([-90 -90],[-20 20],'r',[90 90],[-20 20],'r')

ll=find(UHdir>180);
WSdiff=PSDc1(1:end-1,19)-UHc1(1:end-1,19);[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb(UHdir(ll)-360,WSdiff(ll),5,-180,0);hold on;plot(X,Y,'gd--')
WSdiff=WXTS_TS-UHc1(:,19);[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb([UHdir(ll)-360],WSdiff(ll),5,-180,0);hold on;plot(X,Y,'rd--')
WSdiff=WXTP_TS-UHc1(:,19);[X, Y, M, SD, SE, SDM, SEM, nm] = binave_lb([UHdir(ll)-360],WSdiff(ll),5,-180,0);hold on;plot(X,Y,'bd--')
xlabel('reldir');ylabel('Speed diff (m/s)')
title('5deg-bin avergaed WS difference at 15m')



figure;plot(jd_son,PSDdir,'g',jd_pc,WXTP_TD,jd_pc,WXTS_TD,'r',jd_UH,UHdir,'k');grid;
legend('PSD','SCS port','SCS stbd','UH')

%Will use UH for WS/WD

filesave='C:\data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\cleanedtrueWSD';
eval(['save ' filesave ' jd_UH UH_st UHtruedir'])  
clear filesave


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compares PIR/PSP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PIR1=radmet(:,2);
PIR2=radmet(:,3);
PSP1=radmet(:,4);
PSP2=radmet(:,5);
figure;plot(jd_pc,PIR1,'m',jd_pc,PIR2,'g',jd_pc,WXTS_PIR,'b')

%adjusts to PSP/PIT height
zi=600*ones(size(PSDst));
zref_u=11.5;
zref_t=11.5;
zref_q=11.5;
RHscst=relhum([WXTS_Ta QscsS/0.93-0.34 WXTP_Pa]);
PSDc1=coare30vn_ref([PSDst; NaN],15.9,WXTP_Ta,15,RHscst,15,WXTP_Pa+0.12*(15),WXTP_sst,PSP1,PIR1,WXTP_lat,zi,zref_u,zref_t,zref_q,1);

tt=find(PSP1<1);
figure;plot(jd_pc,WXTP_Ta,jd_pc,Tc1,'m--',jd_pc,Td1,'m',jd_pc,Tc2,'g--',jd_pc,Td2,'g');grid
legend('Tair','Tc1','Td1','Tc2','Td2');
figure;plot(jd_pc(tt),WXTP_Ta(tt),jd_pc(tt),Tc1(tt),'m--',jd_pc(tt),Td1(tt),'m',jd_pc(tt),Tc2(tt),'g--',jd_pc(tt),Td2(tt),'g');grid

figure;
subplot(2,2,1);plot(WXTP_Ta(tt),Tc1(tt),'.',[-5 25],[-5 25],'k--');axis([-5 25 -5 25]);grid
subplot(2,2,2);plot(WXTP_Ta(tt),Td1(tt),'.',[-5 25],[-5 25],'k--');axis([-5 25 -5 25]);grid
subplot(2,2,3);plot(WXTP_Ta(tt),Tc2(tt),'.',[-5 25],[-5 25],'k--');axis([-5 25 -5 25]);grid;
subplot(2,2,4);plot(WXTP_Ta(tt),Td2(tt),'.',[-5 25],[-5 25],'k--');axis([-5 25 -5 25]);grid


figure;plot(jd_pc,PSP1,jd_pc,PSP2,'r',jd_pc,WXTP_PSP,'k')

figure;subplot(2,2,1);plot(WXTP_PSP-15.5,PSP1,'.',[0 800],[0 800],'k--');axis([0 800 0 800]);grid
subplot(2,2,2);plot(WXTP_PSP-15.5,PSP2,'.',[0 800],[0 800],'k--');axis([0 800 0 800]);grid
subplot(2,2,3);plot(PSP1,PSP2,'.',[0 800],[0 800],'k--');axis([0 800 0 800]);grid


figure;subplot(2,2,1);plot(WXTS_PIR,PIR1,'.',[200 450],[200 450],'k--');axis([200 450 200 450]);grid
subplot(2,2,2);plot(WXTS_PIR,PIR2,'.',[200 450],[200 450],'k--');axis([200 450 200 450]);grid
subplot(2,2,3);plot(PIR1,PIR2,'.',[200 450],[200 450],'k--');axis([200 450 200 450]);grid