1. function [s_AB,A_AB,A_BA] = Vinventy(a_fi, a_la, b_fi,b_la,a, e2, eprim2)
2.  
3. b = a*sqrt(1-e2);
4. f = 1-b/a;
5.  
6. tan_U1 = tan(a_fi) * (1-f);
7. tan_U2 = tan(b_fi) * (1-f);
8.  
9. U1 = atan(tan_U1);
10. U2 = atan(tan_U2);
11.  
12. L = b_la - a_la;
13. L_1 = L; %pierwsze przyblizenie
14. etykieta = 0;
15.  
16. i = 1;
17.  
18. while etykieta < 1
19.     if i ~= 1
20.         L_1 = L_2;
21.     end;
22.  
23.     sin_sigma = sqrt( (cos(U2)*sin(L_1))^2 + (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(L_1))^2 );
24.     cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(L_1);
25.  
26.     sigma = atan(sin_sigma/cos_sigma);
27.  
28.     sin_alpha = (cos(U1)*cos(U2)*sin(L_1))/(sin_sigma);
29.     cos_2_alpha = 1 - sin_alpha^2;
30.  
31.     cos_2_sigma_m = cos_sigma - (2*sin(U1)*sin(U2))/(cos_2_alpha);
32.  
33.     C = (f/16) * cos_2_alpha * (4 + f*(4 - 3*cos_2_alpha));
34.  
35.     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)));
36.  
37.     if abs(L_2 - L_1)<(10^(-12)) etykieta = 1;  end;
38.     i = i+1;
39. end;
40.  
41.  
42. u_2 = cos_2_alpha * eprim2;
43. A = 1 + (u_2/16384)*(4096 + u_2*(-768+u_2*(320-175*u_2)));
44. B = (u_2/1024)*(256 + u_2*(-128+u_2*(74-47*u_2)));
45.  
46. 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)));
47.  
48. s_AB = b*A*(sigma - delta_sigma);
49.  
50. A_AB = atan( (cos(U2)*sin(L_2))/(cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(L_2)) );
51. A_BA = atan( (cos(U1)*sin(L_2))/(-sin(U1)*cos(U2) + cos(U1)*sin(U2)*cos(L_2)) );
52. A_BA = A_BA + pi;
53.  
54. end