% routine to compare buoy vs modeled wind

% load data
load /home3/ocn/jmanning/atsee/wind.dat
b=wind;
load /home3/ocn/jmanning/wind/mfiles/modwind.dat
% these are winds m/s for 8, 11, and 5 with forecast hour in 2nd col
m=modwind;


% sort data
[sb,i]=sort(b(:,1));
b=b(i,:);
[sm,i]=sort(m(:,1));
m=m(i,:);

% plot e/w wind data for the three buoys
is=0;
for k=[8 11 5];
  is=is+1;
  subplot(3,2,is*2-1);
  plot(b(find(b(:,6)==k),1),b(find(b(:,6)==k),4),'r','linewidth',2);
  hold on;
  % plot only those forecast <= 24 hours
  plot(m(find(m(:,2)<=24),1),m(find(m(:,2)<=24),2+is*2-1),'g.',...
  'markersize',20);
  % plot points forecasted one whole day prior
  plot(m(find(m(:,2)>24),1),m(find(m(:,2)>24),2+is*2-1),'g*');
  title(['Buoy ',num2str(k)],'fontsize',18);
    axis([349 361 -25 25]);
  if k==5;
    legend_old('Obs U','Mod U');
    yd2calab(98)
  else;
    set(gca,'xticklabels',' ');
  end;
end;

% plot n/s wind data for the three buoys
is=0;
for k=[8 11 5];
  is=is+1;
  subplot(3,2,1+is*2-1);
  plot(b(find(b(:,6)==k),1),b(find(b(:,6)==k),5),'r','linewidth',2);
  hold on;
  % plot only those forecast <= 24 hours
  plot(m(find(m(:,2)<=24),1),m(find(m(:,2)<=24),2+is*2),'g.',...
  'markersize',20);
  % plot points forecasted one whole day prior
  plot(m(find(m(:,2)>24),1),m(find(m(:,2)>24),2+is*2),'g*');
  title(['Buoy ',num2str(k)],'fontsize',18);
    axis([349 361 -25 25]);
  if k==5;
    legend_old('Obs V','Mod V');
    yd2calab(98)
  else;
    set(gca,'xticklabels',' ');
  end;
end;
ibc=0;
for buoy=[11 8 5];
  ibc=ibc+1;
  % calculate error as a function of forecast hour
  for j=6:6:36;
    jc=j/6;
    id=find(m(:,2)==j);
    % for each buoy subtract the observed wind of that time
    eji=0;
    for ej=1:length(id);
      % find buoys at this time
      idb=find(b(:,1)==m(id(ej),1)&b(:,6)==buoy);
      if ~isempty(idb);
        eji=eji+1;
        % for each buoy one at a time
        difu(eji)=b(idb(1),4)-m(id(ej),2+ibc*2-1);
        difv(eji)=b(idb(1),5)-m(id(ej),2+ibc*2);
      end;
    end;
    % at this forecast hour caculate means
    mdifu(jc,ibc)=mean(abs(difu));
    mdifv(jc,ibc)=mean(abs(difv));
    clear difu difv;
  end;
end;

% plot error as a function of forecast hour
figure(2);
plot([6:6:36]', mdifu(:,1),'r',...
       [6:6:36]', mdifu(:,2),'g',...
       [6:6:36]', mdifu(:,3),'b','linewidth',3);
legend_old('buoy 11','buoy 8','buoy 5');
hold on;
plot([6:6:36]', mdifv(:,1),'r--',...
       [6:6:36]', mdifv(:,2),'g--',...
       [6:6:36]', mdifv(:,3),'b--','linewidth',3);
ylabel('Obs-Mod difference (m/s)');
xlabel('Forecast Hour')
title('u solid and v dashed')
set(gcf,'position',[360   80   560   420]);