% plot_adcp_int.m
% =================================================
%       PLOT ADCP: INT OR OA
% =================================================
%
% plot adcp data interpolated on a regular xy grid:
% a) using space interpolation and central date, or
% b) using OA interpolation on each of the days adcp
% data is availablE
% time -- is central time
clear all
addpath /home/valery/b2disk/oa_hershel
addpath /home/valery/bin
addpath /home/valery/b2disk/opnml/FEM
addpath /home/valery/b2disk/opnml/basics
% load current adcp4quoddy velocity file
load adcp4quoddy
% OVERLAY GRID if necessary (grd=1)
           pgrd=1
% integrate over depth 
for i=[1:length(ABS_U)]
u(i)=nanmean(ABS_U(:,i));
v(i)=nanmean(ABS_V(:,i));
end
% 1D filter
%      hh=fspecial('average',[1,13]);
%      uf=filter2(hh,u);vf=filter2(hh,v);
% set the beginning and ending points
      n1=785;n2=length(ABS_U);n2=1450;
% cold eddy C1
      [xx,yy]=meshgrid(-66.6:.05:-65.0,30.5:.05:31.6);
% get central lat, lon
      xcent=mean(xx(:)), ycent=mean(yy(:))
% meshgrid C2
%[xx,yy]=meshgrid(-60.5:.05:-59.0,31.0:.05:31.6);
%mesgrid c1
%[xx,yy]=meshgrid(-61.2:.05:-59.2,30.0:.05:31.0);
      uu=griddata(nav(1,n1:n2),nav(2,n1:n2),u(n1:n2),xx,yy);
      vv=griddata(nav(1,n1:n2),nav(2,n1:n2),v(n1:n2),xx,yy);
      h1=plot(nav(1,n1:n2),nav(2,n1:n2),'b-')
      hold on
      h2=vecplot2(xx(:),yy(:),uu(:),vv(:),.5,'cm/sec')
      set(h2,'color',[.0 .5 .5])
% draw grid
       grid
                  if(pgrd==1)
        load /data/cruises/oc415_1/coord_mtx_05new.txt;
        c=coord_mtx_05new;
        pp=axis
%axis([pp]);
% plot the only numbers inside axis
        ff=find(c(:,3)>pp(1) & c(:,3)<pp(2) & c(:,2)>pp(3) & c(:,2)<pp(4));
      for i=[1:length(ff)]         
        text(c(ff(i),3),c(ff(i),2),['   ',num2str(c(ff(i),1))],'fontsize',7);
      end
      hold on
        plot(c(ff,3),c(ff,2),'sm','markersize',7)
      
  end% BATS site coordinates

%       axis([-66.6  -64.9 32.2 33.6])       
       setaxis
       xlabel('Longitude')
       ylabel('Latitude')
       jj=mean(time(n1:n2))+datenum(2005,1,1);
       title(['Eddy A4: interpolated adcp velocity (centered on  ',datestr(jj,0),' GMT)'])
  return

% now OA adcp
% set days from adcp file
       t=(207:218);
       nav=nav';time=time';u=u';v=v';       
   for k=[1:length(t)]
     for j=[1:length(yy)]
         for i=[1:length(xx)]
[ua(k,j,i),va(k,j,i),err(k,j,i)]=corrUV2(xx(j,i),yy(j,i),t(k),nav(n1:n2,1),nav(n1:n2,2),time(n1:n2),u(n1:n2),v(n1:n2));
         end
     end
     k
   end
       load adcp43_c1_oa
%       figure(2)
       ww=sqrt(uu(:).^2+vv(:).^2);
       ff=find(isnan(ww)==1);
       for k=[1:12]
       clf
       uc=ua(k,:,:);vc=va(k,:,:);
       uc=uc(:);uc(ff)=nan;
       vc=vc(:);vc(ff)=nan;
       gg=find(time > t(k) & time < t(k)+1);
       h5=plot(nav(n1:n2,1),nav(n1:n2,2),'-','color',[.8 .8 .8],'linewidth',2)
       axis([-66.7  -64.1 29.7 32.3]);
       hold on
       h3=plot(nav(gg,1),nav(gg,2),'r-')
       h4=vecplot2(xx(:),yy(:),uc,vc,1.0,'cm/sec')
       set(h4,'color',[.0 .5 .5])
       hold off       
       setaxis
       xlabel('Longitude')
       ylabel('Latitude')
       title(['Eddis C5 & A4: interpolated adcp velocity; day ',num2str(t(k)),'; ', datestr(t(k)+datenum(2004,1,1.5),31)])
       uoa=sqrt(uc(:).^2+vc(:).^2);
       ss(k,1)=sqrt(nanmean((uoa-ww).^2))*100;datestr(t(k)+datenum(2004,1,1.5),31)
       ss(k,2)=nanstd(uoa)*100;
       ss(k,3)=nanstd(ww)*100;
       input('?')
       pp=sprintf('adcp10_c1_%i',t(k));
       pp1=['fig2ppm(''' pp '.ppm'')'];
       eval(pp1);

   end
       figure(3)
       plot(t,ss(:,1),'sk-');
       hold on
       plot(t,ss(:,2),'ob-');
       plot(t,ss(:,3),'*r-');
       hold off
       grid
       xlabel('time, day since 1 Jan, 2004');
       ylabel('rmsq difference, cm/sec');
       title('rmsq difference between interp and OA interp adcp velocity')
       legend('rmsq','oa std','int std',4)
       gtext('t0=10 days')
       
       fig2jpeg('rmsq_10.jpg')
       
       