function s12=distance12(bet1,bet2,c,PV,Geod)

%  s12=distance12(bet1,bet2,c,PV,Geod)computes the distances with given latitudes, 
%  Clairaut constant c ,and PV (passing vertex)  
%  assuming that the points lie on the  reference ellipsoid defined by the input Geod. 
if nargin <5, Geod= defaultellipsoid; end

e=Geod(2);
a=Geod(1);
sbet1=sind(bet1);
cbet1=cosd(bet1);
%[sbet1,cbet1]=SinCosNorm(sbet1,cbet1);
sbet2=sind(bet2);
cbet2=cosd(bet2);
%[sbet2,cbet2]=SinCosNorm(sbet2,cbet2); 

u1=real(sqrt(cbet1.^2-c.^2));
u2=real(sqrt(cbet2.^2-c.^2));
sig12=atan2(sbet2,u2)-(-1).^PV.*atan2(sbet1,u1);   

N1=seriesN1(sbet1,c);
N21=seriesN2(sbet1,c);
N22=seriesN2(sbet2,c);
n=5;
E=e^2.^(1:n); 
N1=N1(:,1:n);
N21=N21(:,1:n);
N22=N22(:,1:n);
N2=N22*E'.*u2-N21*E'.*u1.*(-1).^PV;

s12=a*(PV.*(1+N1*E')*pi+(1+N1*E').*sig12+N2);


end

