Salve a tutti, ho provato a trasformare una funzione trovata su
http://www.movable-type.co.uk/scripts/latlong-vincenty.html (http://www.movable-type.co.uk/scripts/latlong-vincenty.html)
che permette di calcolare la distanza di due punti sulla terra.
In fondo alla pagina è riportato il codice JavaScript che svolge questa funzione
Io ho provato a convertirla in pascal ed il risultato dovrebbe essere:
Function DistVincenty(const Lat1, Lon1, Lat2, Lon2: double): Double;
const a=6378137; // semiasse maggiore della terra in m
b=6356752.314245; // semiasse minore della terra in m
var f, // schiacciamento (a-b)/a
l, // differenza di longitudine
lambda, lambdaP, U1, U2, sinU1, sinU2, cosU1, cosU2, sinLambda, cosLambda,
sinSigma, cosSigma, sigma, sinAlpha, cosSqAlpha, cos2SigmaM, c, uSq,
Aa, Bb, DeltaSigma, s
:double;
IterLimit // limite delle iterazioni
:integer;
begin
//f:=(a-b)/a; // calcolo dello schiacciamento
f:= 1/298.257223563;
l:=degtorad(Lon2-Lon1);
U1:=arctan((1-f)*tan(degtorad(Lat1)));
U2:=arctan((1-f)*tan(degtorad(Lat2)));
sinU1:=sin(U1); cosU1:=cos(U1);
sinU2:=sin(U2); cosU2:=cos(U2);
lambda:=l;
IterLimit:=100;
repeat
sinLambda:=sin(lambda); cosLambda:=cos(lambda);
sinSigma:=sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
(cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
if sinSigma=0 then begin
Result:=0; // punti coincidenti
exit;
end;
cosSigma:=sinU1*sinU2 + cosU1*cosU2*cosLambda;
sigma:=arctan2(sinSigma, cosSigma);
sinAlpha:=cosU1*cosU2*sinLambda/sinSigma;
cosSqAlpha:=1-sinAlpha*sinAlpha;
if cosSqAlpha>0 then
cos2SigmaM:=cosSigma - 2*sinU1*sinU2/cosSqAlpha
else cos2SigmaM:=0; // linea equatoriale
C:= f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
lambdaP:=lambda;
lambda:=L + (1-C) * f * sinAlpha *
(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
Dec(IterLimit);
until ((abs(lambda-lambdaP)>1e-12) and (IterLimit>=0));
if IterLimit=100 then
begin
Result:=NaN;
exit;
end;
uSq:=cosSqAlpha * (a*a - b*b) / (b*b);
Aa:= 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
Bb:= uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
deltaSigma:= Bb*sinSigma*(cos2SigmaM+Bb/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
Bb/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
s:= b*Aa*(sigma-deltaSigma);
Result:=s;
end;
Dove è l'errore?????
Ci stò diventando matto!!! :'( :'(
Datemi una mano per cortesia!!
Salutoni a tutti