Für Koordinaten gibt es grundsätzlich 2 Verfahren. Das erste geht von einer idealen Kugelgestalt der Erde aus. Dies ist sicherlich innerhalb Deutschlands und für deine Belange ausreichend genau. Hierbei geht man meist von einem mittleren Planetenradius von 6371 km aus (die Angaben variieren je nach verwendeter Referenz).
Delphi-Quellcode:
Function CalcSimpleGeoDistance (Lat1, Lon1, Lat2, Lon2, aRadius : Double) : Double;
// Calculates the distance between two positions, based on given
// Latitude and Longitude values on an ideal sqhere with a given
// radius; Dr. Jan Schulz, 2007
Var aPart1 : Double;
aPart2 : Double;
aPart3 : Double;
Begin
Lat1 := Lat1 * Pi / 180;
Lon1 := Lon1 * Pi / 180;
Lat2 := Lat2 * Pi / 180;
Lon2 := Lon2 * Pi / 180;
aPart1 := Cos (Lat1) * Cos (Lon1) * Cos (Lat2) * Cos (Lon2);
aPart2 := Cos (Lat1) * Sin (Lon1) * Cos (Lat2) * Sin (Lon2);
aPart3 := Sin (Lat1) * Sin (Lat2);
Result := ArcCos (aPart1 + aPart2 + aPart3) * aRadius;
end;
Leider ist die Erde aber ein >68-eckiger Klumpen. Für genauere Betrachtungen sollte man deshalb berücksichtigen, daß die Erde durch ihre Rotation zu den Polen hin abgeflacht ist. Besonders wenn man über mehrere 1000 km Distanz mit Koordinaten ausrechnen will. Hier gibt es den Algorithmus von Andoyer, der ihn - um 1949 herum - in einem französischen Ephemeridenkalender veröffentlicht hat. Hab es nicht nachgemessen, aber die Abweichung soll für die Distanz Paris-New York im Bereich der quadrierten Abplattung liegen (<50 m). Natürlich nur in Bezug auf ein Referenzellipsoid. Tatsächlich Weglängen durch Abweichungen im Geoid sind nicht berücksichtigt (Abplattung, Referenzellipsoid und Geoid siehe Wikipedia)...
Delphi-Quellcode:
Function CalcAndoyersGeoDistance (Lat1, Lon1, Lat2, Lon2, RadiusEQ, RadiusPol : Double) : Double;
// Calculates a more precise distance between two positions, based on given
// Latitude and Longitude values for a planet with respect to the flattening
// calculated from the polar and equatorial distance; Dr. Jan Schulz, 2007
Var Flattening : Double;
AndoyerF : Double;
AndoyerG : Double;
AndoyerL : Double;
AndoyerC : Double;
AndoyerS : Double;
AndoyerO : Double;
AndoyerR : Double;
H1 : Double;
H2 : Double;
Begin
Lat1 := Lat1 * Pi / 180;
Lon1 := Lon1 * Pi / 180;
Lat2 := Lat2 * Pi / 180;
Lon2 := Lon2 * Pi / 180;
Flattening := CalcGeometricFlattening (RadiusEQ, RadiusPol);
AndoyerF := (Lat1 + Lat2) / 2;
AndoyerG := (Lat1 - Lat2) / 2;
AndoyerL := (Lon1 - Lon2) / 2;
AndoyerC := Sqr (Cos (AndoyerG)) * Sqr (Cos (AndoyerL)) +
Sqr (Sin (AndoyerF)) * Sqr (Sin (AndoyerL));
AndoyerS := Sqr (Sin (AndoyerG)) * Sqr (Cos (AndoyerL)) +
Sqr (Cos (AndoyerF)) * Sqr (Sin (AndoyerL));
AndoyerO := ArcTan ((Sqrt (AndoyerS / AndoyerC)));
AndoyerR := Sqrt (AndoyerS * AndoyerC) / AndoyerO;
H1 := (3 * AndoyerR - 1) / (2 * AndoyerC);
H2 := (3 * AndoyerR + 1) / (2 * AndoyerS);
Result := 2 * AndoyerO * RadiusEQ * ( 1 +
Flattening * H1 * Sqr (Sin (AndoyerF)) * Sqr (Cos (AndoyerG)) -
Flattening * H2 * Sqr (Cos (AndoyerF)) * Sqr (Sin (AndoyerG)));
end;
Function CalcGeometricFlattening (RadiusEQ, RadiusPol : Double) : Double;
// Calculates the geometric flattening of a planet as relative difference
// between the equatorial and polar radius; Dr. Jan Schulz, 2007
Begin
If RadiusEQ > 0 THen Result := (RadiusEQ - RadiusPol) / RadiusEQ
Else Result := 0;
end;
Bei dieser Methode nimmt man für den Äquatorialradius normalerweise 6378.137 km und für den Polradius 6356.752 km. Je nachdem welches Ellipsoid hierbei zu Grunde gelegt wird, weichen die einzelnen Werte voneinander ab.
Liebe Grüße aus ><)))°> Town
Jan