function [ws, wd, wsave, wdave, ws1, wd1, swath, cnt]=extracthr(speed, dir, lat, lon, land, clat, clon, n, cstart)
%
% [ws wd wsave wdave ws1 wd1 swath cnt]=extracthr(speed, dir, lat, lon, land, clat, clon, N, cstart)
% matlab function to extract averaged winds from BYU high resolution 
% wind fields 
%
% outputs: 
%  (2n by 2n) linear averaged wind: ws,wd
%  (2n by 2n) separately averaged speed and direction: wsave, wdave
%  closest wind ws1, wd1
%  L2B cross-track index value
%  count of averaged vectors
%
% use loadBYUhr to read input arrays
% inputs:
%  selected ambiguity winds: speed, dir
%  winds location and land flag: lat, lon, land 
%  sample point location: clat, clon
%  size of average window (2.5 km pixels): n
%  swath location value

% BYU Sept. 2003

latmin=min(min(lat));
latmax=max(max(lat));
dellat=latmax-latmin;
center_lat=latmin+dellat/2.0;

lonmin=min(min(lon));
lonmax=max(max(lon));
dellon=lonmax-lonmin;
if dellon > 180  % switch from [-180 +180] to [0 360]
  ind=find(lon<0)
  lon(ind)=lon(ind)+360;
  lonmin=min(min(lon));
  lonmax=max(max(lon));
  dellon=lonmax-lonmin;
  if (clon < 0) 
    clon=clon+360;
  end
end

err=(lat-clat).^2+(lon-clon).^2/(cos(pi*center_lat/180))^2;
ind=find(err==min(min(err)));
ind=ind(1);

%[clat lat(ind) clon lon(ind) err(ind)]

if (err(ind) > 0.2)  % miss
  disp(['Missed pass',num2str(err(ind)),' ',num2str(clat),' ',
	num2str(lat(ind)),'  ',num2str(clon),' ',num2str(lon(ind))]);
  ws=0; wd=0; ws1=0; wd1=0;
  return;
end

ws1=speed(ind); % closest vector
wd1=dir(ind);   % closest vector

% compute (2n by 2n) nearby average

[i,j]=ind2sub(size(lat),ind);
[M,N]=size(lat);
i1=max(1,i-n);
i2=min(M,i+n);
j1=max(1,j-n);
j2=min(N,j+n);
[X,Y]=meshgrid([i1:i2],[j1:j2]);
[M1 N1]=size(X);

ind2=sub2ind(size(lat),reshape(X,1,M1*N1),reshape(X,1,M1*N1));

% remove any land or zero estimates
ind=find(land(ind2)==1 & speed(ind2)>0);
ind2=ind2(ind);
cnt=length(ind2)
swath=(i+cstart)/10;  % L2B swath index  

% compute vector mean
u=mean(speed(ind2).*cos(dir(ind2)*pi/180));
v=mean(speed(ind2).*sin(dir(ind2)*pi/180));

ws=u*u+v*v;
wd=0.0;
if (ws > 0.0)
  ws=sqrt(u*u+v*v);
  wd=180*atan2(v,u)/pi;
end

wsave=mean(speed(ind2));  % mean wind speed
wdave=meanangle(dir(ind2));    % mean wind dir

return;

function res=meanangle(dirs)
%
% compute mean of directions (in deg)
%
ang=dirs;
ind=find(ang>180);
ang(ind)=ang(ind)-360;
m1=mean(ang);
e1=mean(180-abs(180-abs(ang-m1)));
ind=find(ang<0);
ang(ind)=ang(ind)+360;
m2=mean(ang);
e2=mean(180-abs(180-abs(ang-m2)));

res=m1;
if e2 < e1
  res=m2;
end
return;




