function [s_AB,A_AB,A_BA] = Vinventy(a_fi, a_la, b_fi,b_la,a, e2, eprim2)

b = a*sqrt(1-e2);
f = 1-b/a;

tan_U1 = tan(a_fi) * (1-f);
tan_U2 = tan(b_fi) * (1-f);

U1 = atan(tan_U1);
U2 = atan(tan_U2);

L = b_la - a_la;
L_1 = L; %pierwsze przyblizenie
etykieta = 0;

i = 1;

while etykieta < 1
    if i ~= 1
        L_1 = L_2;
    end;

    sin_sigma = sqrt( (cos(U2)*sin(L_1))^2 + (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(L_1))^2 );
    cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(L_1);

    sigma = atan(sin_sigma/cos_sigma);

    sin_alpha = (cos(U1)*cos(U2)*sin(L_1))/(sin_sigma);
    cos_2_alpha = 1 - sin_alpha^2;

    cos_2_sigma_m = cos_sigma - (2*sin(U1)*sin(U2))/(cos_2_alpha);

    C = (f/16) * cos_2_alpha * (4 + f*(4 - 3*cos_2_alpha));

    L_2 = L + (1-C)*f*sin_alpha*(sigma + C*sin_sigma*(cos_2_sigma_m + C*cos_sigma*(-1+2*cos_2_sigma_m^2)));

    if abs(L_2 - L_1)<(10^(-12)) etykieta = 1;  end;
    i = i+1;
end;


u_2 = cos_2_alpha * eprim2;
A = 1 + (u_2/16384)*(4096 + u_2*(-768+u_2*(320-175*u_2)));
B = (u_2/1024)*(256 + u_2*(-128+u_2*(74-47*u_2)));

delta_sigma = B*sin_sigma * (cos_2_sigma_m  + 1/4*B*(cos_sigma*(-1+2*cos_2_sigma_m^2) - 1/6*B*cos_2_sigma_m*(-3+4*cos_2_sigma_m^2)));

s_AB = b*A*(sigma - delta_sigma);

A_AB = atan( (cos(U2)*sin(L_2))/(cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(L_2)) );
A_BA = atan( (cos(U1)*sin(L_2))/(-sin(U1)*cos(U2) + cos(U1)*sin(U2)*cos(L_2)) );
A_BA = A_BA + pi;

end