function [wdep,grad,theta]=m_toposlope(lat,lon,dx,dy)
%function [wdep,grad,theta]=m_toposlope(lat,lon,dx,dy)
%determines slope and orientation of bathymetry in a swath
% that is oriented by angle theta in an iterative manner
  %theta is direction of topo gradient (pointing shallow) wrt east
lat=lat(:);	lon=lon(:);
dxy=max(dx,dy);
dlon=2*dxy/(111*cos(max(lat)*pi/180));dlat=2*dxy/111;
mnlat=min(lat)-dlat;mxlat=max(lat)+dlat;
mnlon=min(lon)-dlon;mxlon=max(lon)+dlon;
[p,lons,lats]=m_tbase([mnlon,mxlon,mnlat,mxlat]);
theta=zeros(1,length(lat));grad=theta;wdep=theta;
if dx==dy
  nit=1;
else
  nit=4;
end

for j=1:length(lat);
  yg=(lats-lat(j))*111;
  xg=(lons-lon(j))*cos(lat(j)*pi/180)*111;

  isg=yg<dy & yg>-dy & xg<dx & xg>-dx;
  for i=1:nit
    pg=p(isg)/1000;
    x0=xg(isg);	y0=yg(isg);
    A=[ones(size(x0(:))),x0(:),y0(:)];
    b=pg(:);
    c=A\b;
    wdep(j)=c(1);
    theta(j)=atan2(c(3),c(2));
    if theta(j)<0, theta(j)=theta(j)+pi;end
    grad(j)=sqrt(c(2)^2+c(3)^2);
    tth=tan(theta(j));
    cth=cos(theta(j));
    sth=sin(theta(j));
    isg0=isg;
    isg=yg<-xg/tth+dy/sth & yg>-xg/tth-dy/sth & xg>yg/tth-dx/sth & ...
	xg<yg/tth+dx/sth;
    isg1=find(~(yg<-xg/tth+dy/sth));
    isg2=find(~(yg>-xg/tth-dy/sth));
    isg3=find(~(xg>yg/tth-dx/sth));
    isg4=find(~(xg<yg/tth+dx/sth));
%hold on
%    plot(xg(isg1),yg(isg1),'bx');
%    plot(xg(isg2),yg(isg2),'k+');
%    plot(xg(isg3),yg(isg3),'ro');
%    plot(xg(isg4),yg(isg4),'gd');
    
    
%    keyboard
%    close
  end
end
