function difflon12=dlambda12(bet1,bet2,alp1,PV,e)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
%syms x c 
%Geoid=[6378137 0.0818191908426215]
if nargin < 4, e=0.0818191908426215; end
%PV=(alp1>90);
%c=abs(cosd(bv));
sbet1=sind(bet1);
cbet1=cosd(bet1);
[sbet1,cbet1]=SinCosNorm(sbet1,cbet1);
sbet2=sind(bet2);
cbet2=cosd(bet2);
[sbet2,cbet2]=SinCosNorm(sbet2,cbet2);
salp1=sind(alp1);
calp1=cosd(alp1);
salp0 = salp1 .* cbet1;
c=salp0;
%calp0 = hypot(calp1, salp1 .* sbet1);
ssig1 = sbet1; somg1 = salp1 .* cbet1 .* sbet1;
csig1 = calp1 .* cbet1; comg1 = csig1;
u1=csig1;
[ssig1, csig1] = SinCosNorm(ssig1, csig1);
calp2 =sqrt((calp1.* cbet1).^2+(cbet2-cbet1).* (cbet1 + cbet2))./ cbet2;

ssig2 = sbet2; somg2 =  salp1 .* cbet1 .* sbet2;
csig2 = calp2 .* cbet2;  comg2 = csig2;
u2=csig2;
[ssig2, csig2] = SinCosNorm(ssig2, csig2);
sig12 = atan2(csig1.* ssig2- (-1).^PV.*ssig1.* csig2, ...
                csig1.* csig2+ (-1).^PV.*ssig1.* ssig2);         
omg12 = atan2(max(comg1 .* somg2 - (-1).^PV.*somg1 .* comg2, 0), ...
                comg1 .* comg2 +  (-1).^PV.*somg1 .* somg2);
 sig1=atan2(sbet1,u1);
 sig2=atan2(sbet2,u2);
 
 %testing
  %test
 
 
 
 %testing
dsig1=c./(1-c.^2).*tan(sig1);
dsig2=c./(1-c.^2).*tan(sig2);
dsig12=dsig2- (-1).^PV.*dsig1;

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

dw12=1./(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:2:2*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=diff.*cbet1.*calp1;

end

