function lon12c=longitude12(bet1,bet2,c,PV,e)
%  lon12c=longitude12(bet1,bet2,c,PV,e) computes the longitude difference 
%  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, e=0.0818191908426215; end

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);
omg12=atan2(c.*sbet2,u2)-(-1).^PV.*atan2(c.*sbet1,u1);
           

M1=seriesM1(sbet1,c);
M21=seriesM2(sbet1,c);
M22=seriesM2(sbet2,c);
n=5;
E=e^2.^(1:n); 
M2=M22*E'.*u2-M21*E'.*u1.*(-1).^PV;
lon12=PV.*(1+c.*(M1*E'))*pi+omg12+c.*(M1*E'.*sig12+M2);
lon12(bet1==0 & c==1)=pi*sqrt(1-e^2);
lon12c=rad2deg(lon12);

end

