clear;
          %  L    p   t  e
EOK200 = [ 0.900 740 +15 0 ];
Aga6   = [ 0.550 730 -15 0 ];
Aga70  = [ 0.633 740 -15 0 ];

fyz = [  %   D     p  t   e   dakomer
        14212.836	765	 0.9	32 Aga70
        15359.111	726	 0.0   4 Aga70
         6505.387	737	-0.8	25 Aga6
         7651.662	748	-1.7	47 EOK200
         8798.242	759	-2.6	19 Aga6
         9944.517	720	-3.5	40 EOK200
        11090.793	731	-4.4	12 Aga70
        12237.373	742	-5.3	34 EOK200
        13383.648	753	-6.2   5 Aga70
        14530.229	764	-7.1  27 Aga6
     ];
     
mat = [ %XYZ            A              B             vs    vc    k   phi    azim
      1005758 703610 630.378 1000147 719098 508.823	1.597 1.018	0.20 50.4672 282.448
      1011264 718189 760.393 1012402 714946 108.310	1.662 1.617	0.13 50.4009 101.720
    	1016771 712766 190.365 1004657 710795 407.840	1.726 1.215	0.20 50.3582 1.691
    	1002276 707344 320.380 1016912 706643 707.349	1.791 1.814	0.20 50.4938 169.745
    	1007783 701922 450.374 1009167 702491 306.858	1.856 1.413	0.20 50.4511 194.902
    	1013289 716500 580.367 1001422 718340 606.366	1.921 1.012	0.20 50.3849 343.592
    	1018795 711078 710.382 1013677 714188 205.875	1.986 1.611	0.13 50.3422 321.185
    	1004301 705656 140.355 1005932 710036 505.383	1.051 1.210	0.20 50.4778 242.085
    	1009808 700233 270.348 1018187 705885 104.892	1.116 1.808	0.20 50.4350 206.575
    	1015314 714811 400.363 1010442 701733 404.401	1.181 1.407	0.20 50.3688 61.993
     ];

% ========================================================================= 
R = 6380000;

for ( i=1:size(fyz,1) )
  dd(i,1) = ddBS( fyz(i,5),fyz(i,6),fyz(i,7),fyz(i,8),  fyz(i,2),fyz(i,3),fyz(i,4) );
  if( fyz(i,5:8)==EOK200)
  dd(i,3) = 270.54 - 77978.3 * fyz(i,2) / (273.15+fyz(i,3)) / 740;
  end
  if( fyz(i,5:8)==Aga6)
  dd(i,3) = 309.2 - 109.46 * fyz(i,2) / (273.2+fyz(i,3));
  end
  if( fyz(i,5:8)==Aga70)
  dd(i,3) = 308.6 - 107.9 * fyz(i,2) / (273.2+fyz(i,3));
  end
  D = fyz(i,1) + fyz(i,1)*dd(i,1)*1e-6;
  % refrakce
  k = mat(i,9);
  L = D - (2*k-k^2) * D^3 /24 /R^2;
  % na ref. kouli
  phi = mat(i,10)*pi/180;
  azim = mat(i,11)*pi/180;
  a = 6377397.15508;
  e1 = 0.00667437223;
  e2 = 0.006719218798;
  N = a /sqrt(1-e1*sin(phi)^2);
  Ra = N /(1+e2*cos(phi)^2*cos(azim)^2);
  HA = mat(i,3)+mat(i,7);
  HB = mat(i,6)+mat(i,8);
  l = sqrt(  (L^2 -(HA-HB)^2) /(1+HA/Ra) /(1+HB/Ra)  );
  s = l + l^3 /24 /R^2;
  % na elipsoid
  se = s - HB*D/R *e2 *cos(phi)^2 *cos(azim)^2;
  % do JTSK
  RA = sqrt(mat(i,1)^2+mat(i,2)^2);
  RB = sqrt(mat(i,4)^2+mat(i,5)^2);
  XAB = (mat(i,1)+mat(i,4))/2;
  YAB = (mat(i,2)+mat(i,5))/2;
  RAB = sqrt(XAB^2+YAB^2);
  dd(i,2) = se *( ro2m(RA)+4*ro2m(RAB)+ro2m(RB) ) /6;
end

dd