Формулы для расчета азимута и расстояния

[Леха]

Я новичок в геодезии, поэтому многое здесь для меня непонятно...((
Если можно, то напишите плиз формулы по которым ссылка http://www.garmin.com.ua/tools/calc.php рассчитывает азимут и расстояние между двумя точками.
Если формулы по которым считает сайт неизвестны, то может быть у вас есть свои, уже проверенные? )
Заранее спасибо!

[bingeomap]

http://gis-lab.info/qa/angles-rhumb.html

А вообще та было бы лучше, смотреть на книгу "Высшая геодезия" или "Сфероидическая геодезия".

[Леха]

спасибо, попробую посчитать )

[Franklin1967]

Код по формулам из книжки "Теория фигуры Земли". Считает азимут и дальность в любом направлении на планете с любыми параметрами эллипсоида, включая через полюс

Код: Выделить всё

 function Set3dPoint(X,Y,Z : Double) : T3dPoint;
 begin
  result.X:=X;
  result.Y:=Y;
  result.Z:=Z;
 end;

// ПРЯМАЯ И ОБРАТНАЯ ЗАДАЧИ
// "внешние" переменные заполняемые функцией CalcDelta (см.ниже)
// для анализа в обратной геодезической задаче

function _CalcDelta(Ellps: TEllipsoidProperty;P1:T3dpoint;Dist,Azim:extended;IsCalcK:boolean;var P0: T3dPoint) : T3dPoint;
var  T2, N2, U2, v2, nu,_v, _u, L,
     v_c, vc2, vc3, t4,
     V, cosp1,  t    : extended;
     FCTask : extended;
begin
 FCTask:=sqrt(Ellps.one_es)/Ellps.a;
 t:=tan(P1.X); t2:=sqr(t); t4:=sqr(t2);
 nu:=sqrt(Ellps.es/Ellps.one_es)*cos(P1.X);
 N2:=nu*nu;
 _v:=Dist*sin(Azim); v2:=sqr(_v);
 _u:=Dist*cos(Azim); U2:=sqr(_u);
 cosp1:=1/cos(P1.X);
 V:=sqrt(1+N2);
 v_c:=V*FCTask; vc2:=sqr(v_c); vc3:=power(v_c,3);
 // "внешние" переменные для анализа в обратной геодезической задаче
 // b1
 P0.B := v_c*(1+N2);
 // l1
 P0.L :=v_c*cosp1;
 if iscalcK then exit;

 //L5
 L:=-P0.L*vc3*t*(1+3*t2+N2)/3;
 //dB
 result.X:=_u-0.5*FCTask*(V*sqr(_v)*t+(1+N2)*U2*N2*FCTask *(3*t/FCTask-t2*_u+_u))+
 FCTask*FCTask/6*((1+N2)*v2*(1+3*T2+N2-9*N2*T2)*(0.25*(1+N2)*t*v2*FCTask-_u)+
   P0.B*t*U2*3*(N2*U2*FCTask-(4-13*N2+3*T2*(2-3*N2))*_v/6) +
   P0.B*V*v2*_u*((1+15*t2*(2+3*t2)) *v2/4 - (2+15*t2*(1+t2))*U2)*FCTask/5);
  result.X:= result.X*P0.B+L*N2;
 // dL
 result.Y:=_v*(1+v_c*t*_u)+vc2*((1+3*t2+N2)*U2*_v-t2*V2*_v+v_c*t*(2+3*t2+N2)*U2*_u*_v+
  0.2*vc2*_v*(t2*(1+3*t2)*V2*V2+(2+15*t2+15*t2*t2) *U2*U2-(1+20*t2+30*t2*t2)*V2*U2))/3;
 result.Y:=result.Y*P0.L+L*(V2*_v*_u+N2);
 // dA
 result.Z:=v_c*_v*(t+ 0.5*v_c*(1+2*t2+N2)*_u)+
  vc3*t*((5+6*t2+N2-4*N2*N2) *U2*_v -(1+2*t2+N2)*V2*_v +
   v_c *((1.25+7*t2+6*t4+1.5*N2+2*N2*T2) *U2*_u*_v - (0.25+5*t2+6*t4+0.5*N2+2*N2*T2)*(V2*_v*_u+N2)) +
   vc2*t*((1+20*t2+24*t4)*V2*V2*_v -(58+280*t2+240*t4) *V2*_v*U2 +(61+180*t2+120*t4) *_v*U2*U2)/20)/6;
end;

// ОБРАТНАЯ ЗАДАЧА
// X - dist  Y-A12  z-A21
function ReversGeoTask(Ellps: TEllpsType;const P1,P2 : T3dPoint): t3dPoint;
var DeltaL,DeltaB,MinB,
    MinL,aLen,dA12 : double;
    D,D0,D1,P0 : T3dPoint;
    E : TEllipsoidProperty;
    I : byte;
begin
 E:=Ellipsoids[Ellps];
 D0.X:=P2.X-P1.X;
 D0.Y:=P2.Y-P1.Y;
 MinL:=MaxDouble;
 MinB:=MinL;
 D.X:=D0.X; D.Y:=D0.Y;
 result.Y:=0;
 i:=0; P0:=Set3dPoint(0,0,0);
 dA12:=0; aLen:=0;
 while i<254 do
 begin
   _CalcDelta(E,P1, aLen, dA12, true,P0);
   if D.X=0 then D.X:=D.X+1e-10;
   dA12:=Arctan(P0.b*D.Y/(P0.l*D.X));
   if cos(dA12)=0 then break;
   aLen:=D.X/(P0.b*cos(dA12));
   inc(i);
   D1:=_CalcDelta(E,P1, aLen, dA12,false,P0);
   DeltaB:=Abs(D1.X-D0.X);
   DeltaL:=Abs(D1.Y-D0.Y);
  if (DeltaB<MinB) and (DeltaL<MinL) then
  begin
   MinB:=DeltaB;
   MinL:=DeltaL;
   result.X:=aLen;
   result.Y:=dA12;
  end;
  D.X:=D.X+0.5*DeltaB*(1-2*byte(D.X<D1.X));
  D.Y:=D.Y+0.5*DeltaL*(1-2*byte(D.Y<D1.Y));
  if Abs(D.X)>Abs(5*D0.X) then Break;
 end;
 if result.X<0 then result.Y:=-result.Y;
 result.X:=Abs(result.X);
 if d0.x<=0 then result.y:=pi-result.y;
 if (d0.x>0) and (result.y<0) then  result.y:=2*pi+result.y;
end; 

Параметр эллипсоида указывается в структуре

Код: Выделить всё

TEllipsoidProperty = packed record
  eName   : string;
  a       : double;  // major axis or radius if es=0
  b       : double;  // minor axis or radius if es=0
  e       : double;  // eccentricity
  es      : double;  // e ^ 2
  one_es  : double;  // 1 - e^2
  rone_es : double;  // 1/(1 - e^2)
  ra      : double;  // 1/A
  Rf      : double;  // references
  Code    : string;
 end;

где:
элипсоид WGS84 представляется как
(ename:'WGS 84';a:6378137.0;b:6356752.31424518;e: 0.08181919;es: 0.00669437999014111;
one_es: 0.99330562000985889;rone_es:1.0067394967422762251591434067861;ra:0.15678559428874E-6;
rf: 0.00335281066474748;code:'WGS84'),

Красовский представляется как
es: 0.00672267002233207;one_es: 0.993277329977668;rone_es: 1.00676817019722;ra: 0.156779424519173E-6;rf: 0.00336700336700337;code:'intl'),(ename:'(CK-42) Krassovsky 1942';a:6378245;b:6356863.01877305;e: 0.081813334;