function difflon12=dlongitude12(bet1,bet2,c,PV,e)

%  difflon12=dlongitude12(bet1,bet2,c,PV,e)computes the derivative with respect to
%  Clairaut constant (c) of 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 < 4, 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));    
            
sig1=atan2(sbet1,u1);
sig2=atan2(sbet2,u2);
sig12=sig2-(-1).^PV.*sig1;

dsig12=c./(1-c.^2).*2.*sin(sig12)./(cos(sig2-sig1)+cos(sig2+sig1));
dw12=dsig12./c;

du1=-c./u1;
du2=-c./u2;
M1=seriesM1(sbet1,c);
M21=seriesM2(sbet1,c);
M22=seriesM2(sbet2,c);
dM1=seriesdM1(sbet1,c);
dM21=seriesdM2(sbet1,c);
dM22=seriesdM2(sbet2,c);
n=5;
E=e^2.^(1:n); 
diff=PV.*(dM1*E'*pi)+dw12+dM1*E'.*sig12+c.*(M1*E').*dsig12 ...
    +dM22*E'.*u2-dM21*E'.*u1.*(-1).^PV ...
    +c.*(M22*E').*du2-c.*(M21*E').*du1.*(-1).^PV;
%difflon12=rad2deg(diff.*cbet1.*salp1);
difflon12=rad2deg(diff);

end

