- 01
- 02
- 03
- 04
- 05
- 06
- 07
- 08
- 09
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
function Vincenty(Lat1, Lon1, Lat2, Lon2: Extended): Extended;
const // Параметры эллипсоида:
a = 6378245.0;
f = 1 / 298.3;
b = (1 - f) * a;
EPS = 0.5E-30;
var
APARAM, BPARAM, CPARAM, OMEGA, TanU1, TanU2,
Lambda, LambdaPrev, SinL, CosL, USQR, U1, U2,
SinU1, CosU1, SinU2, CosU2, SinSQSigma, CosSigma,
TanSigma, Sigma, SinAlpha, Cos2SigmaM, DSigma : Extended;
begin
lon1 := lon1 * (PI / 180);
lat1 := lat1 * (PI / 180);
lon2 := lon2 * (PI / 180);
lat2 := lat2 * (PI / 180); //Пересчет значений координат в радианы
TanU1 := (1 - f) * Tan(lat1);
TanU2 := (1 - f) * Tan(lat2);
U1 := ArcTan(TanU1);
U2 := ArcTan(TanU2);
SinCos(U1, SinU1, CosU1);
SinCos(U2, SinU2, CosU2);
OMEGA := lon2 - lon1;
lambda := OMEGA;
repeat //Начало цикла итерации
LambdaPrev:= lambda;
SinCos(lambda, SinL, CosL);
SinSQSigma := (CosU2 * SinL * CosU2 * SinL) +
(CosU1 * SinU2 - SinU1 * CosU2 * CosL) *
(CosU1 * SinU2 - SinU1 * CosU2 * CosL);
CosSigma := SinU1 * SinU2 + CosU1 * CosU2 * CosL;
TanSigma:= Sqrt(SinSQSigma) / CosSigma;
if TanSigma > 0 then Sigma := ArcTan(TanSigma)
else Sigma := ArcTan(TanSigma) + Pi;
if SinSQSigma = 0 then SinAlpha := 0
else SinAlpha := CosU1 * CosU2 * SinL / Sqrt(SinSQSigma);
if (Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))) = 0 then Cos2SigmaM := 0
else Cos2SigmaM:= CosSigma - (2 * SinU1 * SinU2 / (Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))));
CPARAM:= (f / 16) * Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha)) *
(4 + f * (4 - 3 * Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))));
lambda := OMEGA + (1 - CPARAM) * f * SinAlpha * (ArcCos(CosSigma) +
CPARAM * Sin(ArcCos(CosSigma)) * (Cos2SigmaM + CPARAM * CosSigma *
(-1 + 2 * Cos2SigmaM * Cos2SigmaM)));
until Abs(lambda - LambdaPrev) < EPS; // Конец цикла итерации
USQR:= Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha)) *(a * a - b * b) / (b * b);
APARAM := 1 + (USQR / 16384) * (4096 + USQR * (-768 + USQR * (320 - 175 * USQR)));
BPARAM := (USQR / 1024) * (256 + USQR * (-128 + USQR * (74 - 47 * USQR)));
DSigma := BPARAM * SQRT(SinSQSigma) * (Cos2SigmaM + BPARAM / 4 *
(CosSigma * (-1 + 2 * Cos2SigmaM * Cos2SigmaM) - BPARAM / 6 * Cos2SigmaM *
(-3 + 4 * SinSQSigma) * (-3 + 4 * Cos2SigmaM * Cos2SigmaM)));
Result := b * APARAM * (Sigma - DSigma);
end;
{ ©Drkb::04255 }
Алгоритм расчёта километража между двумя точками на земной поверхности методом Винсенти, найден в drkb3.0. Там же весь этот ГК уместился в одной строчке:
distance=sqrt(pow((lon1 - lon2)*111*COS(lat2/57.295781), 2) + pow((lat1) - lat)*111, 2));
, чудноо... :)
Follow us!