function [s12 alp2 sigm]=distance12a(bet1,bet2,alp1,PV,Geod)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
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);

salp1=sind(alp1);
calp1=cosd(alp1);
c = salp1 .* cbet1;
salp2=c./cbet2;
calp2=sqrt(cbet2.^2-c.^2)./cbet2;
alp2=rad2deg(atan2(salp2,calp2));
calp2 =sqrt((calp1.* cbet1).^2+(cbet2-cbet1).* (cbet1 + cbet2))./ cbet2;
calp2a =sqrt((calp1.* cbet1).^2+(sbet1-sbet2).* (sbet1 + sbet2))./ cbet2;
calp2(abs(bet1)<45)=calp2a(abs(bet1)<45);
[salp2, calp2] = SinCosNorm(salp2, calp2);
%salp2=c./cbet2;
%calp2=sqrt(cbet2.^2-c.^2)./cbet2
%u1=cbet1.*calp1
%u2=cbet2.*calp2
u1=cbet1.*calp1;
u2=cbet2.*calp2;
%u2=real(sqrt(cbet2.^2-c.^2))
%u1=cbet1.*cosd(180-90.0003629131842)
%u2=cbet2.*cosd(180-156.392157374217)

ssig1 = sbet1; csig1=u1;%csig1 = calp1 .* cbet1;
[ssig1, csig1] = SinCosNorm(ssig1, csig1);
ssig2 = sbet2; csig2=u2;%csig2 = calp2 .* cbet2;
[ssig2, csig2] = SinCosNorm(ssig2, csig2);
sig12 = atan2(csig1.* ssig2- (-1).^PV.*ssig1.* csig2, ...
                csig1.* csig2+ (-1).^PV.*ssig1.* ssig2);          
            
sig2=atan2(ssig2,csig2);
sig1=atan2(ssig1,csig1);
%sig12 =sig2-(-1).^PV.*sig1;
%
%somg1=c.*sbet1;
%comg1=real(sqrt(cbet1.^2-c.^2));
%w1=atan2(somg1,comg1);
%somg2=c.*sbet2;
%comg2=real(sqrt(cbet2.^2-c.^2));
%w2=atan2(somg2,comg2);
%w12=w2+(-1).^PV*w1
%z1=cbet1*sbet2-sbet1*cbet2*cos(w12)+1i*cbet2*sin(w12)
%sig12
%sig121=phase(cbet1.*cbet2.*cos(w12)+sbet1*sbet2+1i*abs(z1))
%sig121=acos(cbet1.*cbet2.*cos(w12)+sbet1*sbet2)

%[sig12 sig2 sig1]
%sig2=atan2(sbet2,u2);
%sig1=atan2(sbet1,u1);
%sig12=sig2-(-1).^PV.*sig1;

N1=seriesN1(sbet1,c);
N21=seriesN2(sbet1,c);
N22=seriesN2(sbet2,c);
n=5;
E=e.^(2:2:2*n); 
N1=N1(:,1:n);
N21=N21(:,1:n);
N22=N22(:,1:n);
%N2=bet1*0;
%if any(~PV)
%[u2(~PV) u1(~PV) N22(~PV,:)*E'.*u2(~PV)-N21(~PV,:)*E'.*u1(~PV)]
%    N2(~PV)=N22(~PV,:)*E'.*u2(~PV)-N21(~PV,:)*E'.*u1(~PV);
%end
%if any(PV) 
%   N2(PV)=N22(PV,:)*E'.*u2(PV)+N21(PV,:)*E'.*u1(PV)
%end
N2=N22*E'.*u2-N21*E'.*u1.*(-1).^PV;

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

sigm=[sig12 ssig1 csig1 ssig2 csig2];
%S=a*(sigma+N1*E'.*sigma+N2*E'.*u);
%S(vertex)=a*(1+N1(vertex,:)*E')*pi/2.*sign(sbet(vertex));

end

