
function [u,v,u_eq,v_eq]=geovel_from_sla(lon,lat,sla)

% function [u,v,u_eq,v_eq]=geovel_from_sla(lon,lat,sla);
%
% Calculate geostrophic velocity anomaly (m/s)
%       from sea level anomaly (m).
% u, v are f-plane velocities.  u_eq, v_eq are
%       higher-order solutions of a Taylor series
%       expansion for off-equatorial distance
%       (one term more than a beta-plane solution).
% These can be combined as in Lagerloef et al (1999)
%       with a weighting function W=exp(-lat^2/(2.2)^2),
%       u_combined=W*u_eq + (1-W)*u.
%       
% Rick Lumpkin, NOAA/AOML/PhOD

if (length(lon)==size(lon,2))
    lon=lon';
end;
if (length(lat)==size(lat,2))
    lat=lat';
end;

[area,dx,dy]=sw_area(lon,lat);
O=7.2921e-05;
latmat=lat*ones(1,length(lon));
f=2*O*sin(2*pi*latmat/360);
q=find(abs(latmat)<.001);f(q)=NaN;
g=9.81;

[dxE,dyE]=gradient(sla);
dxE=dxE./dx;
dyE=dyE./dy;

u=-g*dyE./f;
v=g*dxE./f;

% equatorial band
rE=6378000;
beta=2*O/rE*cos(2*pi*latmat/360);
% Lagerloef et al 1999:  u=u_0+ lat*u_1 + lat^2*u_2 + ...
%   same for sla, v;   L=1.1132e5 meters, dx eta_0=dy eta_0=0
%   u_n=-g/beta/L*dy eta_(n+1), v_n=g/beta/L*dx eta_(n+1)
L=1.1132e5;  u_eq=NaN*ones(size(u)); v_eq=u_eq;
% polyfit grad eta in band 5S to 5N
qlat=find(abs(lat)<=5);
ordr=3;  % Lagerloef et al: ordr=3
for j=1:length(lon);
    dxe=dxE(qlat,j); dye=dyE(qlat,j);
    utemp=zeros(size(dxe));
    vtemp=utemp;
    if length(find(isfinite(dxe+dye))>ordr)
        q=find(isfinite(dxe+dye));
        lattemp=lat(qlat(q)); 
        dxe=dxe(q); dye=dye(q);
        warning off MATLAB:polyfit:PolyNotUnique
        [px,s]=polyfit(lattemp,dxe,ordr);  % dxe=p(1)*lat^ordr + ... +last(p)
        [py,s]=polyfit(lattemp,dye,ordr);
        warning on MATLAB:polyfit:PolyNotUnique
        % dxe ~ px(1)*lat^ordr + px(2)*lat^(ordr-1) + ... + px(ordr+1)
        %   = Z0 + Z1*theta + Z2*theta^2 + Z3*theta^3
        for i=1:ordr;   % ignore term ordr+1, e.g., mean value of dxe, dye
            utemp=utemp+(lat(qlat).^(ordr-i))     .*( ...
                -g./beta(qlat,j)/L* py(i) );
            vtemp=vtemp+(lat(qlat).^(ordr-i))    .*( ...
                g./beta(qlat,j)/L* px(i) );
        end;
        u_eq(qlat,j)=utemp;
        v_eq(qlat,j)=vtemp;
    end;
end;

if (0)  % plot results
  iin=500:2:700;jin=400:2:600;
  clf;patchmap(lon(jin),lat(iin),sla(iin,jin),[-.2 .2]);
  hold on;
  quiver(lon(jin),lat(iin),4*u(iin,jin),4*v(iin,jin),0);
  hold off;
end;

