% calcmeandrift
%
% Calculate pseudo-Eulerian means of the 
% drifter observations in bins.
% NumObs, NumObsT: drifter days per square degree.

clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (1)  % world
    filein='world_modified';
    Lon=[-180 180]; Lat=[-73 85]; 
    dLon=1; dLat=1;
    fileout='world_modified';
elseif (0)  % entire Atlantic
    filein='atlantic';
    Lon=[-100 25]; Lat=[-70 70];
    dLon=1; dLat=1;     % grid spacing for mean fields
    fileout='atl'; % output file: meandrift_*
elseif (0)  % tropical Atlantic
    filein='trop_atlantic';
    Lon=[-100 20]; Lat=[-25 25];  
    dLon=1; dLat=1;
    fileout='tropatl';
elseif (0)  % tropical Indian
    filein='trop_indian';
    Lon=[30 130]; Lat=[-25 25]; 
    dLon=1; dLat=1;
    fileout='tropind';
elseif (0)  % China
    filein='china';
    Lon=[90 150]; Lat=[0 40]; 
    dLon=1; dLat=1;
    fileout='china';
elseif(0)    % interamerican seas
    filein='atlantic';
    Lon=[-100 -60];Lat=[6 32];
    dLon=.5; dLat=.5;
    fileout='carib';
elseif (0)  % NAC/Mann eddy
    filein='atlantic';
    Lon=[-65 -30]; Lat=[30 55];
    dLon=.5; dLat=.5;
    fileout='nac';
elseif (0)  % test subregion
    filein='test';
    Lon=[-72 -60];Lat=[10 22];
    dLon=1;dLat=1;
    fileout='test';
elseif (0)  % western tropical Pacific
    %filein='west_pac_tropical';
    filein='pacific_93-03';
    Lon=[118 180];Lat=[-30 30];
    dLon=1;dLat=1;
    fileout='wpac';
elseif (0)  % eastern tropical Pacific
    %filein='east_pac_tropical';
    filein='pacific_93-03';
    Lon=[-180 -70];Lat=[-30 30];
    dLon=1;dLat=1;
    fileout='epac';
else    % Hawai
    filein='hawaii';
    Lon=[-180 -130];Lat=[5 35];
    dLon=1; dLat=.5;
    fileout='hawaii';
end;

rotate=1;   % set to 1 to rotate bins in direction
            % of semimajor axis
            
Tlag=5/365.24; % prior Lag.integral time scale for 
            % Gauss-Markov inversion (years)
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if (rotate)
  % define the bin width in the direction of the
  % semimajor and semiminor principal axes
  dX=(dLon+dLat);   % semimajor (degrees)
  dY=(dLon+dLat)/2;   % semiminor (degrees)
  dbin=sqrt(dX.^2+dY.^2);
else
  dbin=dLat; dX=dLon; dY=dLat;
end;

% Within each box, perform a least-squares fit of
% sinusoids onto the individual observations
% to identify the annual and semiannual cycles at
% the same time as the mean values Ubar, Vbar, Tbar.
%
% The annual cycle of SST can be reassembled as:
%  T_annual= T_ann(:,:,1,1)*cos(2*pi*t/365.24) +
%		T_ann(:,:,2,1)*sin(2*pi*t/365.24);
%  error bars are sqrt(sums( (:,:,:,2) components))
%
Lon=min(Lon):dLon:max(Lon);
Lat=min(Lat):dLat:max(Lat);
Ubar=NaN*ones(length(Lat),length(Lon));
Vbar=Ubar;Tbar=Ubar;
eUbar=Ubar; eVbar=Vbar;eTbar=Tbar;
Ubar_Inc=Ubar;Vbar_Inc=Ubar;Tbar_Inc=Ubar;
%dTbar=Ubar; edTbar=Ubar;
UGbar=Ubar;VGbar=Ubar;
eUGbar=Ubar; eVGbar=Vbar;
Ubar_Inc=Ubar;Vbar_Inc=Ubar;Tbar_Inc=Ubar;
NumObs=Ubar;NumObsT=Ubar;
T_ann=NaN*ones(length(Lat),length(Lon),2,2);
U_ann=T_ann;V_ann=T_ann;
T_san=T_ann;U_san=T_ann;V_san=T_ann;
UG_ann=T_ann;VG_ann=T_ann;
UG_san=T_ann;VG_san=T_ann;
%dT_ann=T_ann;dT_san=T_ann;
Season=Ubar;BinAng=Ubar;

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

for i=1:length(Lat);
%for i=71:71;
  fprintf('%d of %d\n',i,length(Lat));
  %eval(['load ',filein,' id lon lat sst u v t u_g v_g d_sst']);
  eval(['load ',filein,' id lon lat sst u v t u_g v_g']);
  % extract data within +-dbin/2 degrees (latitude)
  q=find(lat<Lat(i)-dbin/2 | lat>Lat(i)+dbin/2);
  truncatedrift;
  
  t2=t;t2=t2/365.24;
  t2=t2-floor(t2); % fractional year
  
  for j=1:length(Lon);
    % identify data within +-dbin/2 degrees (longitude)
    if (rotate)
      q=find(lon>=Lon(j)-dbin/2 & lon<=Lon(j)+dbin/2);
      q2=find(isfinite(u(q)+v(q)));
      if (length(q2)>12)
        % calculate angle of semimajor axis
        BinAng(i,j)=.5*atan2( ...
	      4*median( u(q(q2)).*v(q(q2)) ), ...
	      2*(median( u(q(q2)).*u(q(q2)) ) ...
	        -median( v(q(q2)).*v(q(q2)) )));
        %[BinAng(i,j),ex]=pca( ...
        %    u(q(q2))-mean(u(q(q2))), ...
        %    v(q(q2))-mean(v(q(q2))));
        
        % calculate priors for the GM fit
        [y2,u_ap,ex]=fitsine(t(q(q2))/365.24, ...
            u(q(q2)),[1 2],std(u(q(q2))),Tlag);
        [y2,v_ap,ex]=fitsine(t(q(q2))/365.24, ...
            v(q(q2)),[1 2],std(v(q(q2))),Tlag);
      end; % if length(q2)>12
      q2=find(isfinite(u_g(q)+v_g(q)));
      if (length(q2)>12)
        [y2,ug_ap,ex]=fitsine(t(q(q2))/365.24, ...
            u_g(q(q2)),[1 2],std(u_g(q(q2))),Tlag);
        [y2,vg_ap,ex]=fitsine(t(q(q2))/365.24, ...
            v_g(q(q2)),[1 2],std(v_g(q(q2))),Tlag);
      end;  % if length(q2)>12
      q2=find(isfinite(sst(q)));
      if (length(q2)>12)
          [y2,t_ap,ex]=fitsine(t(q(q2))/365.24, ...
              sst(q(q2)),[1 2],std(sst(q(q2))),Tlag);
      end;
      %q2=find(isfinite(d_sst(q)));
      %if (length(q2)>12)
      %    [y2,dt_ap,ex]=fitsine(t(q(q2))/365.24, ...
      %        d_sst(q(q2)),[1 2],std(d_sst(q(q2))),Tlag);
      %end;
      
      % 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 % not rotating
        q=find(lon>=Lon(j)-dbin/2 & lon<=Lon(j)+dbin/2);
        u_ap=[]; v_ap=[]; 
        ug_ap=[]; vg_ap=[];
        t_ap=[]; dt_ap=[];
    end;
    
    q2=find(isfinite(sst(q)));
    NumObsT(i,j)=length(q2)*del_t/dX/dY;
    if (length(q2)>12)
      % fit a mean and annual cycle
      [y,x,ex]=fitsine( ...
          t(q(q2))/365.24,sst(q(q2)),[1 2], ...
          std(sst(q(q2))),Tlag,t_ap);
      
      T_ann(i,j,1,1)=x(1); T_ann(i,j,2,1)=x(3);
      T_ann(i,j,1,2)=ex(1);T_ann(i,j,2,2)=ex(3);
      T_san(i,j,1,1)=x(2); T_san(i,j,2,1)=x(4);
      T_san(i,j,1,2)=ex(2);T_san(i,j,2,2)=ex(4);
      Tbar(i,j)=x(length(x)); eTbar(i,j)=ex(length(ex));
      Tbar_Inc(i,j)=mean(sst(q(q2)));
    end;
    q2=find(isfinite(u(q)+v(q)));
    NumObs(i,j)=length(q2)*del_t/dX/dY;
    if (length(q2)>12)
      % fit a mean and annual cycle
      [y,x,ex]=fitsine( ...
          t(q(q2))/365.24,u(q(q2)),[1 2], ...
          std(u(q(q2))),Tlag,u_ap);
      U_ann(i,j,1,1)=x(1); U_ann(i,j,2,1)=x(3);
      U_ann(i,j,1,2)=ex(1);U_ann(i,j,2,2)=ex(3);
      U_san(i,j,1,1)=x(2); U_san(i,j,2,1)=x(4);
      U_san(i,j,1,2)=ex(2);U_san(i,j,2,2)=ex(4);
      Ubar(i,j)=x(length(x)); eUbar(i,j)=ex(length(ex));
      Ubar_Inc(i,j)=mean(u(q(q2)));
      
      [y,x,ex]=fitsine( ...
          t(q(q2))/365.24,v(q(q2)),[1 2], ...
          std(v(q(q2))),Tlag,v_ap);
      V_ann(i,j,1,1)=x(1); V_ann(i,j,2,1)=x(3);
      V_ann(i,j,1,2)=ex(1);V_ann(i,j,2,2)=ex(3);
      V_san(i,j,1,1)=x(2); V_san(i,j,2,1)=x(4);
      V_san(i,j,1,2)=ex(2);V_san(i,j,2,2)=ex(4);
      Vbar(i,j)=x(length(x)); eVbar(i,j)=ex(length(ex));    
      Vbar_Inc(i,j)=mean(v(q(q2)));
     
      Season(i,j)=mean(exp(2*pi*t2(q(q2))*sqrt(-1)));
    end;    % >12 finite values of u
    
    %q2=find(isfinite(d_sst(q)));
    %if (length(q2)>12)
    %    % fit a mean and annual cycle
    %    [y,x,ex]=fitsine( ...
    %      t(q(q2))/365.24,d_sst(q(q2)),[1 2], ...
    %      std(d_sst(q(q2))),Tlag,dt_ap);
    %        
    %    dT_ann(i,j,1,1)=x(1); dT_ann(i,j,2,1)=x(3);
    %    dT_ann(i,j,1,2)=ex(1);dT_ann(i,j,2,2)=ex(3);
    %    dT_san(i,j,1,1)=x(2); dT_san(i,j,2,1)=x(4);
    %    dT_san(i,j,1,2)=ex(2);dT_san(i,j,2,2)=ex(4);
    %    dTbar(i,j)=x(length(x)); 
    %    edTbar(i,j)=ex(length(ex));
    % end;   % >12 finite values of d_sst
  
     q2=find(isfinite(u_g(q)+v_g(q)));
     if (length(q2)>12)
      
      [y,x,ex]=fitsine( ...
          t(q(q2))/365.24,u_g(q(q2)),[1 2], ...
          std(u_g(q(q2))),Tlag,ug_ap);
      UG_ann(i,j,1,1)=x(1); UG_ann(i,j,2,1)=x(3);
      UG_ann(i,j,1,2)=ex(1);UG_ann(i,j,2,2)=ex(3);
      UG_san(i,j,1,1)=x(2); UG_san(i,j,2,1)=x(4);
      UG_san(i,j,1,2)=ex(2);UG_san(i,j,2,2)=ex(4);
      UGbar(i,j)=x(length(x)); eUGbar(i,j)=ex(length(ex));
      
      [y,x,ex]=fitsine( ...
          t(q(q2))/365.24,v_g(q(q2)),[1 2], ...
          std(v_g(q(q2))),Tlag,vg_ap);
      VG_ann(i,j,1,1)=x(1); VG_ann(i,j,2,1)=x(3);
      VG_ann(i,j,1,2)=ex(1);VG_ann(i,j,2,2)=ex(3);
      VG_san(i,j,1,1)=x(2); VG_san(i,j,2,1)=x(4);
      VG_san(i,j,1,2)=ex(2);VG_san(i,j,2,2)=ex(4);
      VGbar(i,j)=x(length(x)); eVGbar(i,j)=ex(length(ex));
    end; % >12 finite values of u_g
    
  end; % loop over Lon
  patchmap2(Lon,Lat,Tbar);pause(.1);
end;   % loop over Lat
clear i j q q2 y del_t t2 rotate lonny;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate optimally-interpolated values
% on the same Lon, Lat grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% smoothing scale: Detlef's Eulerian length scales
load stammerscales

% allocate memory for interpolated fields
UbarI=NaN*ones(size(Ubar));eUbarI=UbarI;
VbarI=UbarI;eVbarI=UbarI;TbarI=UbarI;eTbarI=UbarI;
UGbarI=UbarI;VGbarI=VbarI;
eUGbarI=UbarI;eVGbarI=VbarI;
UI_ann=NaN*ones(size(U_ann));
VI_ann=UI_ann;UI_san=UI_ann;VI_san=VI_ann;
TI_ann=UI_ann;TI_san=UI_ann;
UGI_ann=UI_ann;VGI_ann=UI_ann;
UGI_san=UI_ann;VGI_san=UI_ann;

% load 1deg Smith and Sandwell topo, to be used
% for coastline detection (no interpolation across land)
load bathyworld
q=find(z2>0);z2(q)=1;
q=find(z2<0);z2(q)=0;

Lonmat=ones(length(Lat),1)*Lon;
Latmat=Lat'*ones(1,length(Lon));
disp('calculating OI fields')
for i=1:length(Lat);
    fprintf('%d of %d\n',i,length(Lat));
    for j=1:length(Lon);
        Le=(1e3*interp2(lon_ds,lat_ds',Le_ds, ...
            Lon(j),Lat(i),'linear'))/100e3;
        D2= (Lonmat-Lon(j)).^2 + (Latmat-Lat(i)).^2;
        W=exp(-(D2 ./ Le.^2));
        
        % identify where W>.01 but the straight-line path
        % from bin i,j to this bin extends over land.
        % if bin i,j is on land, don't block unless the
        % path is over water first (e.g. interpolation
        % allowed along a coastline, with bin centers on land).
        q=find(W>.01 & W<1);
        for ii=1:length(q);
            x=[Lonmat(i,j),Lonmat(q(ii))];
            y=[Latmat(i,j),Latmat(q(ii))];
            if (abs(diff(x))>abs(diff(y)))
                x2=x(1):.1*diff(x):x(2);
                y2=interp1(x,y,x2);
            else
                y2=y(1):.1*diff(y):y(2);
                x2=interp1(y,x,y2);
            end;
            land=interp2(lon2,lat2,z2,x2,y2);
            if (max(diff(land))>0)  % point blocked by land
                W(q(ii))=0;
            end;
        end;
        
        q=find(W>.01 & isfinite(Ubar+eUbar));
        if (length(q))
            [UbarI(i,j),eUbarI(i,j)]=interp2_opt( ...
                Lonmat(q),Latmat(q),Ubar(q), ...
                Lon(j),Lat(i),Le,eUbar(q));
        end;
        q=find(W>.01 & isfinite(Vbar+eVbar));
        if (length(q))
            [VbarI(i,j),eVbarI(i,j)]=interp2_opt( ...
                Lonmat(q),Latmat(q),Vbar(q), ...
                Lon(j),Lat(i),Le,eVbar(q));
        end;
        q=find(W>.01 & isfinite(Tbar+eTbar));
        if (length(q))
            [TbarI(i,j),eTbarI(i,j)]=interp2_opt( ...
                Lonmat(q),Latmat(q),Tbar(q), ...
                Lon(j),Lat(i),Le,eTbar(q));
        end;
        
        q=find(W>.01 & isfinite(UGbar+eUGbar));
        if (length(q))
            [UGbarI(i,j),eUGbarI(i,j)]=interp2_opt( ...
                Lonmat(q),Latmat(q),UGbar(q), ...
                Lon(j),Lat(i),Le,eUGbar(q));
        end;
        q=find(W>.01 & isfinite(VGbar+eVGbar));
        if (length(q))
            [VGbarI(i,j),eVGbarI(i,j)]=interp2_opt( ...
                Lonmat(q),Latmat(q),VGbar(q), ...
                Lon(j),Lat(i),Le,eVGbar(q));
        end;
        for comp=1:2;
            ua=U_ann(:,:,comp,1);eua=U_ann(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                UI_ann(i,j,comp,1)=ua;
                UI_ann(i,j,comp,2)=eua;
            end;

            ua=V_ann(:,:,comp,1);eua=V_ann(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                VI_ann(i,j,comp,1)=ua;
                VI_ann(i,j,comp,2)=eua;
            end;
            
            ua=T_ann(:,:,comp,1);eua=T_ann(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                TI_ann(i,j,comp,1)=ua;
                TI_ann(i,j,comp,2)=eua;
            end;
            
            ua=U_san(:,:,comp,1);eua=U_san(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                UI_san(i,j,comp,1)=ua;
                UI_san(i,j,comp,2)=eua;
            end;

            ua=V_san(:,:,comp,1);eua=V_san(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                VI_san(i,j,comp,1)=ua;
                VI_san(i,j,comp,2)=eua;
            end;
            
            ua=T_san(:,:,comp,1);eua=T_san(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                TI_san(i,j,comp,1)=ua;
                TI_san(i,j,comp,2)=eua;
            end;
 
            ua=UG_ann(:,:,comp,1);eua=UG_ann(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                UGI_ann(i,j,comp,1)=ua;
                UGI_ann(i,j,comp,2)=eua;
            end;

            ua=VG_ann(:,:,comp,1);eua=VG_ann(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                VGI_ann(i,j,comp,1)=ua;
                VGI_ann(i,j,comp,2)=eua;
            end;
            
            ua=UG_san(:,:,comp,1);eua=UG_san(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                UGI_san(i,j,comp,1)=ua;
                UGI_san(i,j,comp,2)=eua;
            end;

            ua=VG_san(:,:,comp,1);eua=VG_san(:,:,comp,2);
            q=find(W>.01 & isfinite(ua+eua));
            if (length(q))
                [ua,eua]=interp2_opt( ...
                    Lonmat(q),Latmat(q),ua(q), ...
                    Lon(j),Lat(i),Le,eua(q));
                VGI_san(i,j,comp,1)=ua;
                VGI_san(i,j,comp,2)=eua;
            end;
        end;    % loop over components 1,2
        
    end;    % loop over Lon (j)
end;        % loop over Lat (i)
clear Lonmat Latmat *_ds W i j D2 ua eua 
clear Le x ex comp;

% calculate mean Ekman velocities
%[Uek_soc,Vek_soc]=calcmeanekman(Lon,Lat,1);
if (0)
    [Uek,Vek]=calcmeanekman( ...
        [min(Lon)-2,Lon,max(Lon)+2], ...
        [min(Lat)-2,Lat,max(Lat)+2],2);
    Uek=Uek(2:length(Lat)+1,2:length(Lon)+1,:);
    Vek=Vek(2:length(Lat)+1,2:length(Lon)+1,:);
    Uek=mean(Uek,3);Vek=mean(Vek,3);
end;

fprintf('saving file meandrift_%s.mat',fileout);
eval(['save meandrift_',fileout, ...
	' Lon Lat *bar* NumObs* *_ann* *_san*', ...
	' Season dX dY BinAng Tlag']);
