function  wd=w_depth(xi,yi)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  wd=W_DEPTH(xi,yi) computes the water depth at the point with

%  longitude xi and latitude yi using the basic 3.5 km bathymetric 

%  database of D. Greenberg (February 1, 1996, GLOBEC Homepage).

%

%    INPUT:  xi = longitude (decimal degrees, west is negative)

%            yi = latitude  (decimal degrees)

%

%    OUTPUT: wd = water depth (m)

%

%  Assumes xi,yi are vectors. wd=999.9 for those input points 

%  which are located outside of the grid. 

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 11/11/96

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%load gomdepth

load w_depthsf


% set default depth, determine size

wd=999.9.*ones(size(xi));

M=length(xi);



% constants

L=3.75; % approximate length of box

k=pi./180;

r=6370; % earth radius

rk=r.*k;

theta=34.1; % mean counterclockwise rotation of bathymetric database

CR=cos(k.*theta);

SR=cos(k.*theta);



% initiate loop 

for m=1:M;

  

xin=xi(m);

yin=yi(m);



% compute separation distance of input point from all database points

C=cos(k.*yin);

C2=C.^2;

dist=rk.*sqrt((ygom-yin).^2+C2.*(xgom-xin).^2);



% locate closest gridpoint in database array

[Dist,I]=sort(dist);

I1=I(1);

xc=xgom(I1);

yc=ygom(I1);



% compute dxr,dyr using Pin as origin

dx=rk.*C.*(xin-xc);

dy=rk.*(yin-yc);

dxr=CR.*dx+SR.*dy;

dyr=-SR.*dx+CR.*dy;



% skip calculation if Pin outside grid

distmin=L/sqrt(2)+0.1;

if abs(dxr)<distmin & abs(dyr)<distmin

   

% identify which quadrant contains Pin

if dxr>0 & dyr>0

   quad=1;

 elseif dxr<0 & dyr>0

   quad=2;

 elseif dxr>0 & dyr<0

   quad=3;

 elseif dxr<0 & dyr<0

   quad=4;

end



% compute indices for each quadrant

dlat=L.*CR./rk;dlon=L.*SR./(rk.*C);

if quad==1

   J1=I1;J2=J1-1;

   x3=xc-dlon;y3=yc+dlat;

   dist=rk.*sqrt((ygom-y3).^2+C2.*(xgom-x3).^2);

   [Dist,I]=sort(dist);

   J3=I(1);J4=J3-1;

elseif quad==2

   J2=I1;J1=J2+1;

   x4=xc-dlon;y4=yc+dlat;

   dist=rk.*sqrt((ygom-y4).^2+C2.*(xgom-x4).^2);

   [Dist,I]=sort(dist);

   J4=I(1);J3=J4+1;

elseif quad==3

   J3=I1;J4=J3-1;

   x1=xc+dlon;y1=yc-dlat;

   dist=rk.*sqrt((ygom-y1).^2+C2.*(xgom-x1).^2);

   [Dist,I]=sort(dist);

   J1=I(1);J2=J1-1;

elseif quad==4

   J4=I1;J3=J4+1;

   x2=xc+dlon;y2=yc-dlat;

   dist=rk.*sqrt((ygom-y2).^2+C2.*(xgom-x2).^2);

   [Dist,I]=sort(dist);

   J2=I(1);J1=J2+1;

end



% check on logic

J=[J1 J2 J3 J4];

xp=xgom(J);

yp=ygom(J);

dp=dgom(J);

% plot(xp,yp,'+',xin,yin,'o',xp(1),yp(1),'x')

% pause



% compute dx's,dy's from Pin

dx=rk.*C.*(xin-xp);

dy=rk.*(yin-yp);



% set up array and solve for depth at input point

A=zeros(4,4);

for mm=1:4

  A(mm,:)=[1 dx(mm) dy(mm) dx(mm).*dy(mm)];

end

depth=A\dp;

wd(m)=depth(1);



end



end

 