Obliczanie odległości na podstawie długości i szerokości geogr.

Witam,

mam działający kod w Delphi, który oblicza odległość w km na podstawie szerokości i długości geograficznej danych punktów:

function TfrmMain.GetDistanceBetween(long1, lat1, long2, lat2 : Double) : Double;
var
  F,G,L : Double;
  SF, SG, SL,
  CF, CG, CL : Double;
  W1, W2 : Double;
  S, C : Double;
  O,R,D : Double;
  H1, H2 : Double;
  ff : Double;
begin
  ff := 1 / 298.257;
  F := (lat1 + lat2) / 2;
  G := (lat1 - lat2) / 2;
  L := (long1 - long2) / 2;

  SF := Sin(F*Pi/180);
  SG := Sin(G*Pi/180);
  SL := Sin(L*Pi/180);
  CF := Cos(F*Pi/180);
  CG := Cos(G*Pi/180);
  CL := Cos(L*Pi/180);

  W1 := sqr(SG * CL);
  W2 := sqr(CF * SL);
  S := W1 + W2;

  W1 := sqr(CG * CL);
  W2 := sqr(SF * SL);
  C := W1 + W2;
  O := ArcTan(Sqrt(S/C));
  R := Sqrt(S*C) / O;
  D := 2 * O * 6378.14;

  H1 := (3*R-1) / (2*C);
  H2 := (3*R+1) / (2*S);

  W1 := sqr(SF * CG) * H1 * ff + 1;
  W2 := sqr(CF * SG) * H2 * ff;

  result := D * (W1 - W2) * 1.609344
end;

Ale problem jest, kiedy muszę przenieść tą funkcję do procedure w bazie Firebird.

Kod wygłada tak:

CREATE PROCEDURE OBLICZ_ODLEGOSC(
  LONG1 DOUBLE PRECISION,
  LAT1 DOUBLE PRECISION,
  LONG2 DOUBLE PRECISION,
  LAT2 DOUBLE PRECISION)
RETURNS(
  DISTANCE NUMERIC(12, 1))
AS
DECLARE VARIABLE F DOUBLE PRECISION;
DECLARE VARIABLE G DOUBLE PRECISION;
DECLARE VARIABLE L DOUBLE PRECISION;
DECLARE VARIABLE SF DOUBLE PRECISION;
DECLARE VARIABLE SG DOUBLE PRECISION;
DECLARE VARIABLE SL DOUBLE PRECISION;
DECLARE VARIABLE CF DOUBLE PRECISION;
DECLARE VARIABLE CG DOUBLE PRECISION;
DECLARE VARIABLE CL DOUBLE PRECISION;
DECLARE VARIABLE W1 DOUBLE PRECISION;
DECLARE VARIABLE W2 DOUBLE PRECISION;
DECLARE VARIABLE S DOUBLE PRECISION;
DECLARE VARIABLE C DOUBLE PRECISION;
DECLARE VARIABLE O DOUBLE PRECISION;
DECLARE VARIABLE R DOUBLE PRECISION;
DECLARE VARIABLE D DOUBLE PRECISION;
DECLARE VARIABLE H1 DOUBLE PRECISION;
DECLARE VARIABLE H2 DOUBLE PRECISION;
DECLARE VARIABLE ff DOUBLE PRECISION;
begin
  ff = 1 / 298.257;
  F = (lat1 + lat2) / 2;
  G = (lat1 - lat2) / 2;
  L = (long1 - long2) / 2; 
  
  SF = Sin(F * Pi() / 180);
  SG = Sin(G * Pi() / 180);
  SL = Sin(L * Pi() / 180);
  CF = Cos(F * Pi() / 180);
  CG = Cos(G * Pi() / 180);
  CL = Cos(L * Pi() / 180);


  W1 = sqr(SG * CL); -- nie ma SQR
  W2 = sqr(CF * SL); -- nie ma SQR
  S = W1 + W2;

  W1 = sqr(CG * CL); -- nie ma SQR
  W2 = sqr(SF * SL); -- nie ma SQR
  C = W1 + W2;
  O = ATan(Sqrt(S/C));
  R = Sqrt(S*C) / O;
  D = 2 * O * 6378.14;

  H1 = (3*R-1) / (2*C);
  H2 = (3*R+1) / (2*S);

  W1 = sqr(SF * CG) * H1 * ff + 1; -- nie ma SQR
  W2 = sqr(CF * SG) * H2 * ff; -- nie ma SQR

  DISTANCE = D * (W1 - W2) * 1.609344;
  SUSPEND;
END;

W firebirdzie jest funkcja SQRT, ale nie ma SQR. Pytanie do mózgów matematycznych, w jaki sposób mogę sprawić, aby powyższy kod zadziałał?

Sqr(A*B) = (A*B)*(A*B)

Czyli jak dobrze rozumiem:

W1 = (CG * CL)*(CG * CL);

Ja tam z matmy orłem nie jestem, ale zdaje mi się, że SQR(A) = A * A...

Tak, sqr jest to podniesienie liczby do kwadratu (do potęgi drugiej), a potęgowanie to mnożenie liczby przez siebie odpowiednią ilość razy.

Wielkie dzięki :) Działa!

Ja liczę znacznie prostszym wzorem:

function Dist(Lat1,Lon1,Lat2,Lon2:Double):Double;
const dz=12756.274;//średnica Ziemi na równiku [km]
var a,b:Double;
begin
  a:=(Lon2-Lon1)*Cos(Lat1*pi/180);
  b:=(Lat2-Lat1);
  Result:=Sqrt(a*a+b*b)*pi*dz/360;//[km]
end;