
% calcmeaneddy
%
% calculate the mean eddy kinetic energy, time scales, etc.
% from the mean-removed drifter currents.

% Gain: the slope of a best-fit line between
%  u_gp_tot and u_av.  <1 means that u_av is generally
%   weaker than u_gp_tot.
clear all;
filename='world'; % meandrift_* extension

eval(['load meandrift_',filename,' Lon Lat dX dY BinAng']);
dLon=median(diff(Lon));
dLat=median(diff(Lat));
if (~isempty(find(isfinite(BinAng), 1)))
    rotate=1;
    dbin=sqrt(dX.^2+dY.^2);
else
    rotate=0;
    dbin=dLat;
end;

eval(['load eddyterms_',filename,' t']);
del_t=median(diff(t(1:1000)));
clear t;

Upbar=NaN*ones(length(Lat),length(Lon));
Vpbar=Upbar; Tpbar=Upbar;
UGpbar=Upbar; VGpbar=Upbar;
Tlubar=Upbar; Tlvbar=Upbar;
Tdubar=Upbar; Tdvbar=Upbar;
Tlubar_std=Upbar; Tlvbar_std=Upbar;
Tltbar=Upbar; Tdtbar=Upbar;
NumObs=Upbar; NumObsT=Upbar; NumObsLag=Upbar;
UpVpbar=Upbar; UpTpbar=Upbar; VpTpbar=Upbar;
UGpVGpbar=Upbar;


for i=1:length(Lat); 
  fprintf('%d of %d\n',i,length(Lat));
  tic
  eval(['load eddyterms_',filename,' lon lat t id']);
  eval(['load eddyterms_',filename,' up vp sstp Tl* Td*']);
  eval(['load eddyterms_',filename,' *_gp_tot *_gp *_av']);
  q=find(lat<Lat(i)-dbin/2 | lat>Lat(i)+dbin/2);
  truncatedrift;
  
  for j=1:length(Lon); 
    if (rotate)
      % identify data within rotated rectangle
      x=(lon-Lon(j)).*cos(BinAng(i,j))+ ...
          (lat-Lat(i)).*sin(BinAng(i,j));
      y=-(lon-Lon(j)).*sin(BinAng(i,j))+ ...
          (lat-Lat(i)).*cos(BinAng(i,j));
      q=find(x>=-dX/2 & x<=dX/2);
      q2=find(y(q)>=-dY/2 & y(q)<=dY/2);
      q=q(q2);
    else
      q=find(lon>=Lon(j)-dLon/2 & lon<=Lon(j)+dLon/2);
    end;

    Tlubar(i,j)=nanmean(Tlu(q));
    Tlvbar(i,j)=nanmean(Tlv(q));
    Tdubar(i,j)=nanmean(Tdu(q));
    Tdvbar(i,j)=nanmean(Tdv(q));
    NumObsLag(i,j)=length(isfinite(Tlu(q))) ...
	  *del_t/dX/dY;
    
    if (length(q)>12)
      Tlubar_std(i,j)=nanstd(Tlu(q));
      Tlvbar_std(i,j)=nanstd(Tlv(q));
      NumObs(i,j)=length(isfinite(up(q)))*del_t/dX/dY;
      Upbar(i,j)=sqrt(nanmean(up(q).^2));
      Vpbar(i,j)=sqrt(nanmean(vp(q).^2));
      UGpbar(i,j)=sqrt(nanmean(u_gp(q).^2));
      VGpbar(i,j)=sqrt(nanmean(v_gp(q).^2));
      UpVpbar(i,j)=nanmean(up(q).*vp(q));
      UGpVGpbar(i,j)=nanmean(u_gp(q).*v_gp(q));
      q2=find(isfinite(sstp(q)));
      if (length(q2)>12);
        NumObsT(i,j)=length(q2)*del_t/dX/dY;
        Tpbar(i,j)=sqrt(nanmean(sstp(q).^2));
        UpTpbar(i,j)=nanmean(up(q).*sstp(q));
        VpTpbar(i,j)=nanmean(vp(q).*sstp(q));
      end;
      %q2=find(isfinite(u_av(q)+u_gp_tot(q)));
      %if (length(q2)>12);
          %c=corrcoef(u_av(q(q2)),u_gp_tot(q(q2)));
          %Corr_u(i,j)=c(1,2);
          %c=corrcoef(v_av(q(q2)),v_gp_tot(q(q2)));
          %Corr_v(i,j)=c(1,2);
          %c=polyfit(u_gp_tot(q(q2)),u_av(q(q2)),1);
          %Gain_u(i,j)=c(1);
          %c=polyfit(v_gp_tot(q(q2)),v_av(q(q2)),1);
          %Gain_v(i,j)=c(1);
          
      %end;
    end;
  end; % loop over Lon (index j)
  
  
  pcolor(Lon,Lat,sqrt(Upbar.^2+Vpbar.^2));
  caxis([0 .5]);shading flat;
  pause(.01);
  toc
end;   % loop over Lat (index i)
clear i j q q2 y rotate;

% calculate number of independent observations
% in each bin, using an integral timescale of
% 2 days where isnan(Tlubar+Tlvbar) and NumObs>0.
tlu=Tlubar+Tlvbar;
q=find(isnan(tlu) & isfinite(NumObs));
tlu(q)=2;
Nind=NumObs*dX*dY ./ ...
        weightedsmooth(tlu,NumObsLag);

% temp: using Lag time for velocity here ...
%    longer timescale for temperature?
NTind=NumObsT*dX*dY ./ ...
        weightedsmooth(tlu,NumObsLag);
clear tlu q

fprintf('saving file meaneddy_%s.mat\n',filename);
eval(['save meaneddy_',filename, ...
        ' Lon Lat *bar *bar_std NumObs* N*ind']);

    
    
% synthda - synthesize drifters and AVISO altimetry.
 

eval(['load eddyterms_',filename, ...
        ' lon lat *_gp_* *_av t id']);
q=find(isnan(u_gp_tot+u_av));truncatedrift; 
eval(['load meandrift_',filename,' UGbarI VGbarI;'])
LonB=Lon; LatB=Lat;

dLon=median(diff(Lon));
dLat=median(diff(Lat));
C=NaN*ones(length(Lat),length(Lon));
A=C;Uc=C;Vc=C;
NumPts=C;

for i=1:length(Lat)
    fprintf('%d of %d\n',i,length(Lat));
    q=find(lat>=Lat(i)-dLat/2 & lat<=Lat(i)+dLat/2);
    for j=1:length(Lon)
        q2=find(lon(q)>=Lon(j)-dLon/2 & lon(q)<=Lon(j)+dLon/2);
        qq=q(q2);
        
        if (~isempty(qq))

        NumPts(i,j)=length(qq);
        % calculate mean dot product {{v_dg cdot vs}}
        vdv=mean(dot([u_gp_tot(qq) v_gp_tot(qq)]', ...
            [u_av(qq) v_av(qq)]'));
        
        % calculate correlation coefficient
        C(i,j)=vdv./( sqrt( ...
            mean( u_gp_tot(qq).^2 + v_gp_tot(qq).^2 ) * ...
            mean( u_av(qq).^2 + v_av(qq).^2) ));
        
        % calculate mean dot product {{vs cdot vs}}
        vsdvs=mean( u_av(qq).^2+v_av(qq).^2 );
        % calculate amplitude adjustment A
        A(i,j)=( vdv-dot([mean(u_gp_tot(qq)) mean(v_gp_tot(qq))], ...
            [mean(u_av(qq)) mean(v_av(qq))]) ) ./ ...
            (vsdvs - dot([mean(u_av(qq)) mean(v_av(qq))], ...
            [mean(u_av(qq)) mean(v_av(qq))]) );
        % calculate time-mean geostrophic velocity v_c
        Uc(i,j)=interp2(LonB,LatB',UGbarI,Lon(j),Lat(i),'linear') ...
            -A(i,j)*mean(u_av(qq));
        Vc(i,j)=interp2(LonB,LatB',VGbarI,Lon(j),Lat(i),'linear') ...
            -A(i,j)*mean(v_av(qq));
        
        end; % if (length(qq))
    end;    % loop over Lon (j)
    if (length(isfinite(A))>2)
        patchmap2(Lon,Lat,A);pause(.1);
    end;
end;    % loop over Lat (i)

P=quiver(Lon,Lat,5*Uc,5*Vc,0);set(P,'color','k');
axis([min(Lon) max(Lon) min(Lat) max(Lat)])
overlay;asprat;

eval(['load meaneddy_',filename]);
eval(['save meaneddy_',filename, ...
        ' Lon Lat *bar *bar_std NumObs* N*ind A C Uc Vc NumPts']);

