%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set path to include my .m files
%
%path(path,'../matlab_cen');
%
% Load observational data
%
clear all
OPATH='../OBS/';
UVOF=[OPATH,'UV_INTEG.m3d'];

[UVO,gridname,year,ncol,gmt]=read_obs(UVOF);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% separate out moorings and adcp
% UVOB, UVOE, UVOI = data for moorings B, E, I respectively (south to north)
% UVO_adcp = adcp data
Aind = find( (round(UVO(:,2)) == 999999) & (round(UVO(:,3)) == 999999) );
Bind = find( (round(UVO(:,2)) == 49280) & (round(UVO(:,3)) == 92859) );
Eind = find( (round(UVO(:,2)) == 137162) & (round(UVO(:,3)) == 152816) );
Iind = find( (round(UVO(:,2)) == 239828) & (round(UVO(:,3)) == 198418) );

% 2003 (less decimal places?)
%Bind = find( (round(UVO(:,2)) == 49444) & (round(UVO(:,3)) == 92859) );
%Eind = find( (round(UVO(:,2)) == 137490) & (round(UVO(:,3)) == 153270) );
%Iind = find( (round(UVO(:,2)) == 239910) & (round(UVO(:,3)) == 197960) );

UVOA = UVO(Aind,:);
UVOB = UVO(Bind,:);
UVOE = UVO(Eind,:);
UVOI = UVO(Iind,:);

UVO_temp = UVO;
UVO_temp([Bind' Eind' Iind'],:) = NaN;
ADCPind = find(isfinite(UVO_temp(:,1))); 
UVO_adcp = UVO_temp(ADCPind,:);
clear UVO_temp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load model results
%
for isim=1:3
   MPATH=['../CRUISE_Q5.',num2str(isim),'/'];
   UVMF=[MPATH,'UV_INTEG-Q5.m3d'];
   [UVEM(isim,:,:),gridname,year,ncol,gmt]=read_obs_trunc(UVMF);
end
gridname='/home/cstock/Data/meshes/stw2/stw2';
gridstruct=loadgrid(gridname);

% filter out any points that are outside the domain by finding the closest 
% model data comparison for each observation point.  If this different is
% larger than "eps" (set to 2m here) then remove for purposes of comparison

UVO_temp = UVO;
ncomp = size(UVO,1);
count = 0;
eps = 2; 
for nn = 1:ncomp;
  aa = find( (abs(UVEM(1,:,2) - UVO(nn,2)) < eps) & (abs(UVEM(1,:,3) - UVO(nn,3)) < eps));
  if isempty(aa);
    count = count+1;
    UVO_temp(nn,:) = [NaN NaN NaN NaN NaN NaN NaN NaN];
  end
end

ind = find(isfinite(UVO_temp(:,1)));
UVOf = UVO_temp(ind,:);

% separate the moorings and the ADCP for model misfit data - we need
% this for the stellwagon mesh because some of the data lies outside
% the grid boundary. UVOf_adcp is the filtered version of the ADCP.

Aind_m = find( (round(UVEM(1,:,2)) == 999999) & (round(UVEM(1,:,3)) == 99999) );
Bind_m = find( (round(UVEM(1,:,2)) == 49280) & (round(UVEM(1,:,3)) == 92859) );
Eind_m = find( (round(UVEM(1,:,2)) == 137162) & (round(UVEM(1,:,3)) == 152816) );
Iind_m = find( (round(UVEM(1,:,2)) == 239828) & (round(UVEM(1,:,3)) == 198418) );

% 2003 (OC391)
%Bind_m = find( (round(UVEM(1,:,2)) == 49444) & (round(UVEM(1,:,3)) == 92859) );
%Eind_m = find( (round(UVEM(1,:,2)) == 137490) & (round(UVEM(1,:,3)) == 153270) );
%Iind_m = find( (round(UVEM(1,:,2)) == 239910) & (round(UVEM(1,:,3)) == 197960) );

for isim = 1:3
  UVEMA(isim,:,:) = UVEM(isim,Aind_m,:);
  UVEMB(isim,:,:) = UVEM(isim,Bind_m,:);
  UVEME(isim,:,:) = UVEM(isim,Eind_m,:);
  UVEMI(isim,:,:) = UVEM(isim,Iind_m,:);
end

UVEM_temp = UVEM;
UVEM_temp(:,[Bind_m Eind_m Iind_m],:) = NaN;
ADCPind_m = find(isfinite(UVEM_temp(1,:,1)));
UVEM_adcp = UVEM_temp(:,ADCPind_m,:);
clear UVEM_temp;

% adcp observation with 'out of bounds' points filtered
UVOf_adcp = UVOf(ADCPind_m,:);

% Plot time-series of observational and model results
%
%tallfigure;
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do U,V ADCP comparisons if you have data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare ADCP_u
if ~isempty(UVOf_adcp)
figure(1)
clf
for isim=1:3
   iip = (isim-1)*2 + 1;
   subplot(3,2,iip);
   hold on
   hp=plot(UVOf_adcp(:,1),UVOf_adcp(:,4),'b-');
   nmod=length(UVEM_adcp(isim,:,1));
% Take this out, there are some points that are ouside of the domain
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOf_adcp(:,4)-reshape(UVEM_adcp(isim,:,4),nmod,1);
   VM=UVOf_adcp(:,5)-reshape(UVEM_adcp(isim,:,5),nmod,1);
   hp=plot(UVEM_adcp(isim,:,1),UM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs U','Mod U');
   ylabel('Velocity (m/s)');
   title(['ADCP u ',num2str(isim),': rms error=',num2str(rms(UVEM_adcp(isim,:,4)+sqrt(-1)*UVEM_adcp(isim,:,5))),' m/s']);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end

% compare ADCP_v
   subplot(3,2,iip+1) 
   hold on 
   hp=plot(UVOf_adcp(:,1),UVOf_adcp(:,5),'b-','MarkerSize',6);
   nmod=length(UVEM_adcp(isim,:,1));
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOf_adcp(:,4)-reshape(UVEM_adcp(isim,:,4),nmod,1);
   VM=UVOf_adcp(:,5)-reshape(UVEM_adcp(isim,:,5),nmod,1);
   hp=plot(UVEM_adcp(isim,:,1),UM,'r--','MarkerSize',6);
   set(hp,'LineWidth',1.0)
   legend('Obs V','Mod V');
   ylabel('Velocity (m/s)');

   nobs = size(UVOf_adcp,1) 
   title(['v, bias:', num2str(sum(UVEM_adcp(isim,:,4)+i*UVEM_adcp(isim,:,5))/nobs)]);

   if isim == 3
     xlabel('Time (ydgmt0)');
   end
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc adcp.eps

end
% End ADCP comparison

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do mooring comparisons if you have data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mooring A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(Aind)

figure(2)
clf
for isim=1:3
   iip = (isim-1)*2 + 1;
   subplot(3,2,iip);
   hold on
   hp=plot(UVOB(:,1),UVOB(:,4),'b-');
   nmod=length(UVEMA(isim,:,1));
% Take this out, there are some points that are ouside of the domain
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOA(:,4)-reshape(UVEMB(isim,:,4),nmod,1);
   VM=UVOA(:,5)-reshape(UVEMB(isim,:,5),nmod,1);
   hp=plot(UVEMA(isim,:,1),UM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs U','Mod U');
   ylabel('Velocity (m/s)');
   title(['Mooring A, u',num2str(isim),': rms error=',num2str(rms(UVEMA(isim,:,4)+sqrt(-1)*UVEMA(isim,:,5))),' m/s']);
   if (isim == 3)
     xlabel('Time (ydgmt0)');
   end
   subplot(3,2,iip+1)
   hold on
   hp=plot(UVOA(:,1),UVOA(:,5),'b-');
   hp=plot(UVEMB(isim,:,1),VM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs V','Mod V');
   ylabel('Velocity (m/s)');
   title('v');
   if (isim == 3)
     xlabel('Time (ydgmt0)');
   end
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc mooring_a.eps

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mooring B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(Bind)
figure(3)
clf
for isim=1:3
   iip = (isim-1)*2 + 1;
   subplot(3,2,iip);
   hold on
   hp=plot(UVOB(:,1),UVOB(:,4),'b-');
   nmod=length(UVEMB(isim,:,1));
% Take this out, there are some points that are ouside of the domain
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOB(:,4)-reshape(UVEMB(isim,:,4),nmod,1);
   VM=UVOB(:,5)-reshape(UVEMB(isim,:,5),nmod,1);
   hp=plot(UVEMB(isim,:,1),UM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs U','Mod U');
   ylabel('Velocity (m/s)');
   title(['Mooring B, u ',num2str(isim),': rms error=',num2str(rms(UVEMB(isim,:,4)+sqrt(-1)*UVEMB(isim,:,5))),' m/s']);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end

   subplot(3,2,iip+1)
   hold on
   hp=plot(UVOB(:,1),UVOB(:,5),'b-');
   hp=plot(UVEMB(isim,:,1),VM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs V','Mod V');
   ylabel('Velocity (m/s)');
   nobs = size(UVOB,1);
   title(['v, bias:', num2str(sum(UVEMB(isim,:,4)+i*UVEMB(isim,:,5))/nobs)]);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc mooring_b.eps

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mooring E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(Eind)
figure(4)
clf
for isim=1:3
   iip = (isim-1)*2 + 1;
   subplot(3,2,iip);
   hold on
   hp=plot(UVOE(:,1),UVOE(:,4),'b-');
   nmod=length(UVEME(isim,:,1));
% Take this out, there are some points that are ouside of the domain
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOE(:,4)-reshape(UVEME(isim,:,4),nmod,1);
   VM=UVOE(:,5)-reshape(UVEME(isim,:,5),nmod,1);
   hp=plot(UVEME(isim,:,1),UM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs U','Mod U');
   ylabel('Velocity (m/s)');
   title(['Mooring E, u ',num2str(isim),': rms error=',num2str(rms(UVEME(isim,:,4)+sqrt(-1)*UVEME(isim,:,5))),' m/s']);
   if (isim == 3) 
     xlabel('Time (ydgmt0)');
   end

   subplot(3,2,iip+1)
   hold on
   hp=plot(UVOE(:,1),UVOE(:,5),'b-');
   hp=plot(UVEME(isim,:,1),VM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs V','Mod V');
   ylabel('Velocity (m/s)');
   nobs = size(UVOE,1);
   title(['v, bias:', num2str(sum(UVEME(isim,:,4)+i*UVEME(isim,:,5))/nobs)]);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc mooring_e.eps

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mooring I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(Iind)
figure(5)
clf
for isim=1:3
   iip = (isim-1)*2 + 1;
   subplot(3,2,iip);
   hold on
   hp=plot(UVOI(:,1),UVOI(:,4),'b-');
   nmod=length(UVEMI(isim,:,1));
% Take this out, there are some points that are ouside of the domain
%   for i=1:length(UVEM(isim,:,1))
%      ii(i)=find(abs(UVO(:,1)-UVEM(isim,i,1))<0.0001);
%   end
   UM=UVOI(:,4)-reshape(UVEMI(isim,:,4),nmod,1);
   VM=UVOI(:,5)-reshape(UVEMI(isim,:,5),nmod,1);
   hp=plot(UVEMI(isim,:,1),UM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs U','Mod U');
   ylabel('Velocity (m/s)');
   title(['Mooring I, u',num2str(isim),': rms error=',num2str(rms(UVEMI(isim,:,4)+sqrt(-1)*UVEMI(isim,:,5))),' m/s']);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end

   subplot(3,2,iip+1)
   hold on
   hp=plot(UVOI(:,1),UVOI(:,5),'b-');
   hp=plot(UVEMI(isim,:,1),VM,'r--');
   set(hp,'LineWidth',1.0);
   legend('Obs V','Mod V');
   ylabel('Velocity (m/s)');
   nobs = size(UVOI,1);
   title(['v, bias:', num2str(sum(UVEMI(isim,:,4)+i*UVEMI(isim,:,5))/nobs)]);
   if isim == 3
     xlabel('Time (ydgmt0)');
   end
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc mooring_i.eps

end
% end moorings
%
% Vector plots of observational and model results
%
%tallfigure;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% all data within domain plotted as a function of space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(6)
subplot(2,2,1)
plotbnd(gridstruct);
axis([0 5e5 -1e5 4e5]);
axis('square');set(gca,'XTick',[]);set(gca,'YTick',[]);
hold on
% V = UVO(ii,:);
V=UVOf(:,:);
hv=vecplot2(V(:,2),V(:,3),V(:,4),V(:,5),0.5,'m/s',230000.,-160000.);
set(hv,'color','blue');
title(verbatim(['All observations: RMS = ',num2str(rms(V(:,4)+sqrt(-1.0)*V(:,5))),' m/s']));
for isim=1:3
   subplot(2,2,isim+1)
   V=reshape(UVEM(isim,:,:),length(UVEM(isim,:,1)),ncol);
   plotbnd(gridstruct);
   axis([0 5e5 -1e5 4e5]);
   axis('square');set(gca,'XTick',[]);set(gca,'YTick',[]);
   hold on
   V=reshape(UVEM(isim,:,:),length(UVEM(isim,:,1)),ncol);
   hv=vecplot2(V(:,2),V(:,3),V(:,4),V(:,5),0.5,'m/s',230000.,-160000.);
   set(hv,'color','red');
   title(['Quoddy Simulation ',num2str(isim),': rms error=',num2str(rms(V(:,4)+sqrt(-1)*V(:,5))),' m/s']);
end

set(gcf,'PaperPosition',[0.5 0.25 8 10]);
print -depsc all_residuals.eps


