function [range,A12,A21]=distsphr(lat,long,spheroid);
% DIST    Computes distance and bearing between points on the earth
%         using various reference spheroids.
%
%         [RANGE,AF,AR]=distsphr(LAT,LONG) computes the ranges RANGE between
%         points specified in the LAT and LONG vectors (degrees). Forward
%         and reverse bearings (degrees) are returned in AF, AR.
%
%         [..]=DIST(...,'ellipsoid') uses the specified ellipsoid
%         to get distances and bearing. Available ellipsoids are:
%
%         'clarke66'  Clarke 1866
%         'iau73'     IAU 1973
%         'wgs84'     WGS 1984
%         'sphere'    Sphere of radius 6371.0 km
%
%          The default is 'wgs84'.


%Notes: RP (WHOI) 3/Dec/91
%         Mostly copied from BDC "dist.f" routine (copied from ....?), but
%         then wildly modified to bring it in line with Matlab vectorization.

%C GIVEN THE LATITUDES AND LONGITUDES (IN DEG.) IT ASSUMES THE IAU SPHERO
%C DEFINED IN THE NOTES ON PAGE 523 OF THE EXPLANATORY SUPPLEMENT TO THE
%C AMERICAN EPHEMERIS.
%C
%C THIS PROGRAM COMPUTES THE DISTANCE ALONG THE NORMAL
%C SECTION (IN M.) OF A SPECIFIED REFERENCE SPHEROID GIVEN
%C THE GEODETIC LATITUDES AND LONGITUDES OF THE END POINTS
%C  *** IN DEC%C
%C  IT USES ROBBIN'S FORMULA, AS GIVEN BY BOMFORD, GEODESY,
%C FOURTH EDITION, P. 122.  CORRECT TO ONE PART IN 10**8
%C AT 1600 KM.  ERRORS OF 20 M AT 5000 KM.
%C
%C   CHECK:  SMITHSONIAN METEOROLOGICAL TABLES, PP. 483 AND 484,
%C GIVES LENGTHS OF ONE DEGREE OF LATITUDE AND LONGITUDE
%C AS A FUNCTION OF LATITUDE. (SO DOES THE EPHEMERIS ABOVE)
%C
%C PETER WORCESTER, AS TOLD TO BRUCE CORNUELLE...1983 MAY 27
%C

if (nargin<3), spheroid='wgs84'; end;

if (spheroid(1:3)=='sph'),
      A = 6371000.0;
      B = A;
      E = sqrt(A*A-B*B)/A;
      EPS= E*E/(1-E*E);
elseif (spheroid(1:3)=='cla'), 
      A = 6378206.4E0;
      B = 6356583.8E0;
      E= sqrt(A*A-B*B)/A;
      EPS = E*E/(1.-E*E);
elseif(spheroid(1:3)=='iau'),
      A = 6378160.e0;
      B = 6356774.516E0;
      E = sqrt(A*A-B*B)/A; 
      EPS = E*E/(1.-E*E);
elseif(spheroid(1:3)=='wgs'),

%c on 9/11/88, Peter Worcester gave me the constants for the 
%c WGS84 spheroid, and he gave A (semi-major axis), F = (A-B)/A
%c (flattening) (where B is the semi-minor axis), and E is the
%c eccentricity, E = ( (A**2 - B**2)**.5 )/ A
%c the numbers from peter are: A=6378137.; 1/F = 298.257223563
%c E = 0.081819191
      A = 6378137.;
      E = 0.081819191;
      EPS= E*E/(1.-E*E);

else
   error('dist: Unknown spheroid specified!');
end;


NN=max(size(lat));
if (NN ~= max(size(long))), 
   error('dist: Lat, Long vectors of different sizes!');
end

lat=lat*pi/180;      % convert to radians
long=long*pi/180;

PHI1=lat(1:NN-1);    % endpoints of each segment
XLAM1=long(1:NN-1);
PHI2=lat(2:NN);
XLAM2=long(2:NN);

                    % wiggle lines of constant lat to prevent numerical probs.
if (any(PHI1==PHI2)), 
   for ii=1:NN-1,
      if (PHI1(ii)==PHI2(ii)), PHI2(ii)=PHI2(ii)+ 1e-14; end;
   end;
end;
                     % wiggle lines of constant long to prevent numerical probs.
if (any(XLAM1==XLAM2)), 
   for ii=1:NN-1,
      if (XLAM1(ii)==XLAM2(ii)), XLAM2(ii)=XLAM2(ii)+ 1e-14; end;
   end;
end;



%C  COMPUTE THE RADIUS OF CURVATURE IN THE PRIME VERTICAL FOR
%C EACH POINT

xnu=A./sqrt(1.0-(E*sin(lat)).^2);
xnu1=xnu(1:NN-1);
xnu2=xnu(2:NN);

%C*** COMPUTE THE AZIMUTHS.  A12 (A21) IS THE AZIMUTH AT POINT 1 (2)
%C OF THE NORMAL SECTION CONTAINING THE POINT 2 (1)

TPSI2=(1.-E*E)*tan(PHI2) + E*E*xnu1.*sin(PHI1)./(xnu2.*cos(PHI2));
PSI2=atan(TPSI2);

%C*** SOME FORM OF ANGLE DIFFERENCE COMPUTED HERE??

DPHI2=PHI2-PSI2;
DLAM=XLAM2-XLAM1;
CTA12=(cos(PHI1).*TPSI2 - sin(PHI1).*cos(DLAM))./sin(DLAM);
A12=atan((1.)./CTA12);
CTA21P=(sin(PSI2).*cos(DLAM) - cos(PSI2).*tan(PHI1))./sin(DLAM);
A21P=atan((1.)./CTA21P);

%C    GET THE QUADRANT RIGHT

A12 =A12  + pi*( (abs(A21P-A12)<=(pi/2)) & ( (PHI2-PHI1)<0. ) );
A21P=A21P + pi*( (abs(A21P-A12)<=(pi/2)) & ( (PHI2-PHI1)>=0.) );

SSIG=sin(DLAM).*cos(PSI2)./sin(A12);
SIG=asin(SSIG);
SSIGC=-sin(DLAM).*cos(PHI1)./sin(A21P);
SIGC=asin(SSIGC);
A21 = A21P - DPHI2.*sin(A21P).*tan(SIG/2.0);

%C   COMPUTE RANGE

G2=EPS*(sin(PHI1)).^2;
G=sqrt(G2);
H2=EPS*(cos(PHI1).*cos(A12)).^2;
H=sqrt(H2);
TERM1=-SIG.*SIG.*H2.*(1.0-H2)/6.0;
TERM2=(SIG.^3).*G.*H.*(1.0-2.0*H2)/8.0;
TERM3=(SIG.^4).*(H2.*(4.0-7.0*H2)-3.0*G2.*(1.0-7.0*H2))/120.0;
TERM4=-(SIG.^5).*G.*H/48.0;

range=xnu1.*SIG.*(1.0+TERM1+TERM2+TERM3+TERM4);


%C*** CONVERT TO DECIMAL DEGREES

A12=A12*180/pi;
A21=A21*180/pi;





 














