function dist = gcdist_sr(A,reflat,reflon)
% by T. Rossby
% Routine computes great circle distance z in kilometers on a spheroid
% with flattening =1/295=.00339 using the Andoyer-Lambert approximation.
% dist is the distance (km) from a (lat,long) to b (lat,long) in degrees. 
% 6378.2 is the equatorial radius of the earth
% Reference point given by E. Gottlieb for Oleander : 74 W, 40.75 N 
% NOAA's reference point for Route MB given in XBT manual is 41 23' N  75 W
% Glenn Strout gave me a reference for the TSG data as 40.5N 74.0W
% if nargin > 3
%  reflat = 40.75; reflon = -74;
% end
% reflat = 40.75; reflon = -74;

%%-- Create reference point vector (Rossby Lab standard is 40.75 N 74 W)
B = ones(size(A)); B(:,1) = B(:,1)*(reflat); B(:,2)=B(:,2)*(reflon);

rad = 180/3.14159;
lat2 = A(:,1)/rad;
lon2 = A(:,2)/rad;
lat1 = B(:,1)/rad;
lon1 = B(:,2)/rad;
snst = sin(lat1);
snsc = sin(lat2);
cfz = sin(lat1).*sin(lat2)+cos(lat1).*cos(lat2).*cos(lon1-lon2);
z = 1 - cfz.*cfz;
sfz = sqrt(z);
z = atan2(sfz,cfz);
a = snst - snsc;
b = snst + snsc;
delz = .84752e-3*((z+3*sfz)./(1-cfz).*a.*a + (z-3*sfz)./(1+cfz).*b.*b);
dist = 6378.2*(z - delz);
return
