%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;startup;
%
% Load observational data
%
OPATH='../OBS/'
UVOF=[OPATH,'DRIFTERS.m3d']
%
[UVO,gridname,year,ncol,gmt]=read_obs(UVOF);
no=length(UVO(:,1));
ndr=6;ndts=no/ndr;
TO=(reshape(UVO(:,1),ndr,ndts))';
XO=(reshape(UVO(:,2),ndr,ndts))';
YO=(reshape(UVO(:,3),ndr,ndts))';
UO=(reshape(UVO(:,5),ndr,ndts))';
VO=(reshape(UVO(:,6),ndr,ndts))';
%
% Load model results
%
isim=input('Enter the QUODDY4 simulation number (i.e. 1,2,3): ');
MPATH=['../CRUISE_Q4.',num2str(isim),'/'];
UVMF=[MPATH,'DRIFTERS.Q4.m3d']
[UVM,gridname,year,ncol,gmt]=read_obs_trunc(UVMF);
UM=(reshape(UVM(:,5),ndr,ndts))';
VM=(reshape(UVM(:,6),ndr,ndts))';
%
% Remove drifter 6
ndr=5
TO=TO(:,1:ndr);
XO=XO(:,1:ndr);
YO=YO(:,1:ndr);
UO=UO(:,1:ndr);
VO=VO(:,1:ndr);
UM=UM(:,1:ndr);
VM=VM(:,1:ndr);
%
% Truncate for stationary drifters
X=XO;Y=YO;
DX(2:ndts,:)=X(2:ndts,:)-X(1:ndts-1,:);DX(1,:)=DX(2,:);
DY(2:ndts,:)=Y(2:ndts,:)-Y(1:ndts-1,:);DY(1,:)=DY(2,:);
issta=find(DX==0 & DY==0);
TO(issta)=NaN;
XO(issta)=NaN;
YO(issta)=NaN;
UO(issta)=NaN;
VO(issta)=NaN;
UM(issta)=NaN;
VM(issta)=NaN;
%
%  Generate statistics
VECO(1,:,:)=UO';VECO(2,:,:)=VO';VECO(3,:,:)=0;
VECM(1,:,:)=UM';VECM(2,:,:)=VM';VECM(3,:,:)=0;
OCROSSM=cross(VECO,VECM);
OCROSSM=(reshape(OCROSSM(3,:),ndr,ndts))';
SPDO=sqrt(UO.^2+VO.^2);
SPDM=sqrt(UM.^2+VM.^2);
THETA=180.0/acos(-1.0)*asin(OCROSSM./SPDM./SPDO);
thetamean=(nanmean(THETA'))';
%
UE=UO-UM;
VE=VO-VM;
spderr=sqrt(UE.^2+VE.^2);
spderrmean=(nanmean(spderr'))';
gridname='../MESH/bank150';
[in,x,y,z,bnd]=loadgrid(gridname);
%
% Time-domain plots
%
tallfigure;hold on;figone=gcf;
clrli=['rx- ';'b-. ';'ro--';'r^-.';'b-- ';'b-  '];
for i=1:5
   string(i,:)=['Drifter ',num2str(i)];
   subplot(2,1,1);hold on;
   hp=plot(TO(:,i),100.0*spderr(:,i),clrli(i,:));
   subplot(2,1,2);hold on;
   hp=plot(TO(:,i),THETA(:,i),clrli(i,:));
   title(['Drifter ',num2str(i),]);
end
subplot(2,1,1)
axis([163.5 168 0 50.0]);
hp=plot(TO(:,4),100.0*spderrmean,'k-');
set(hp,'LineWidth',2.0)
string(6,:)=['mean err ']
legend(string);
title(verbatim([MPATH,': Speed differences along drifter tracks (cm/s)']));
xlabel('time (yearday)');
ylabel('|\Delta velocity|');
text(164.0,45.0,['mean(|Verr|) = ',num2str(0.1*round(10.0*nanmean(100.0*spderrmean))),' cm/s']);
%
subplot(2,1,2)
axis([163.5 168 -100.0 100.0]);
hp=plot(TO(:,4),thetamean,'k-');
set(hp,'LineWidth',2.0)
string(6,:)=['mean ang ']
legend(string);
title(verbatim([MPATH,': Angle differences along drifter tracks (deg)']));
xlabel('time (yearday)');
ylabel('|\Delta angle|');
text(164.0,80.0,['mean(|\Theta err|) = ',num2str(0.1*round(10.0*nanmean(thetamean))),' deg']);
%
eval(['print -depsc DRIFTERS_CRUISE_Q4.',num2str(isim),'a.ps']);
eval(['!snatch "Figure No. ',num2str(figone),' " DRIFTERS_CRUISE_Q4.',num2str(isim),'a']);
pause;
%
% Vector plots of observational and model results
%
tallfigure;hold on;figtwo=gcf;
%
plotbnd(x,y,bnd);
axis([170000 220000 -220000 -165000]);
axis('square');hold on;
hv=vecplot2(XO,YO,UO,VO,2.0,'cm/s',210000.,-168000.);
set(hv,'color','blue');
title(verbatim(['Observations: RMS = ',num2str(0.1*round(10.0*100.0*rms(UVO(:,5)+sqrt(-1.0)*UVO(:,6)))),'cm/s']))
hv=vecplot2(XO,YO,UM,VM,2.0,'cm/s',210000.,-173000.);
set(hv,'color','red')
xlabel(['Quoddy Simulation ',num2str(isim),': RMS =',num2str(0.1*round(10.0*100.0*rms(UVM(:,5)+sqrt(-1)*UVM(:,6)))),'cm/s']);
%
TI=[164;165;166;167]
for i=1:5
[smin,nclosest]=closestnode(TO(:,i),TO(:,i),TI,TI)
hp=plot(XO(nclosest,i),YO(nclosest,i),'ko');
set(hp,'LineWidth',2.0);
hp=text(XO(nclosest,i),YO(nclosest,i),num2str(TI));
set(hp,'FontWeight','bold');
end
%
eval(['print -depsc DRIFTERS_CRUISE_Q4.',num2str(isim),'b.ps']);
eval(['!snatch "Figure No. ',num2str(figone),' " DRIFTERS_CRUISE_Q4.',num2str(isim),'b']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
