function lon12c=lambda12(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);
calp2a =sqrt((calp1.* cbet1).^2+(cbet2-cbet1).* (cbet1 + cbet2))./ cbet2;
calp2b =sqrt((calp1.* cbet1).^2+(sbet1-sbet2).* (sbet1 + sbet2))./ cbet2;
calp2=(
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(comg1 .* somg2 - (-1).^PV.*somg1 .* comg2, ...
                comg1 .* comg2 +  (-1).^PV.*somg1 .* somg2);
            
           
 %test
 %c=salp1.* cbet1;
 %u1=cbet1.*calp1;
 %u2=real(sqrt(cbet2.^2-c.^2));
 %sig1=atan2(sbet1,u1);
 %sig2=atan2(sbet2,u2);
 %sig12=sig2-(-1).^PV.*sig1;
 %somg1=c.*sbet1;
 %comg1=real(sqrt(cbet1.^2-c.^2));
 %somg2=c.*sbet2;
 %comg2=real(sqrt(cbet2.^2-c.^2));
 
 %omg1=atan2(somg1,comg1)
 %omg2=atan2(somg2,comg2)
 %omg12=omg2-(-1).^PV.*omg1
 %omg12=phase(cos(sig2)+1i*c.*sin(sig2))-phase(cos(sig1)+1i*c.*sin(sig1));
 
 %

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

lon12c=rad2deg(lon12);

end

