지구에서 어느 두 지점의 경도(longitude)와 위도(latitude) 값을 알면 그 사이의 거리를 구할 수 있다?
말은 쉬워 보입니다. 하지만 여러분들은 이 공식을 유도해본 적 있나요? 그리 쉽지만은 않을 겁니다. 그러나 간단한 수학적 접근만 하면 충분히 유도할 수 있습니다. 이 공식을 구하기 위해선 좌표계 변환(coordinate transform)과 삼각함수(Trigonometrical function)를 이용하면 됩니다. 여기서는 직접 다음 소개될 공식을 유도하지 않습니다. 유도는 여러분께 맡기고요. 저는 이 유도된 공식을 쓸 수 있는 C와 C# 프로그램만 소개할까 합니다.
이해를 돕기위해 그림을 하나 그렸습니다. (무척 힘들었다. 무려 20분 걸렸음 ㅡㅡ;;)
자, 위의 그림은 지구입니다. 노란표시된 부분 위쪽이 우리나라이고 밑에는 호주입니다.(믿거나 말거나) 그럼 우리는 우리나라의 한 좌표의 경도, 위도값(lon1, lat1)과 호주의 한 좌표의 경도, 위도값(lon2, lat2)를 안다고 가정했을때 구면에서의 두 좌표의 거리 d를 구하는게 목표입니다.
잠깐 넘겨집고 간다면 경도(longitude)는 동경 0도~180도, 서경 0~180도 까지 범위를 가집니다. 동경일때는 +로 하고 서경일때는 -라고 생각합시다. 또 위도는 북위 0~90도, 남위 0~90도 까지가 범위입니다. 마찬가지로 북위는 +, 남위는 -라고 합시다. 이렇게 가정한 하에 프로그램을 제작할 것입니다.
자 그럼 공식을 소개하죠.
자 생각보다 어렵지 않습니다. 유도하는데 시간이 걸리겠지만요. 여기서 한가지 가정은 지구가 정확한 구가 아니라는 겁니다. 그러나 잘 알다시피 지구는 표면도 거칠 뿐 아니라 정확한 구가 아닙니다. 적도 쪽으로 약간 더 찌그러진 타원체이지요. 아주 정확한 거리를 계산하려면 이 방법을 쓰면 안됩니다. 타원체일때를 고려할려면 다른 공식을 계산해서 써야겠지요. 그건 나중에 소개하고요.(언제? ㅋ)
그럼 저 공식을 이용해서 만든 C#코드를 공개하겠습니다. 참고로 이 코드는 Codeproject.com에서 주어왔어요. ^^
C# 코드 (Language : cpp)
using System.Text;
public class CDistanceBetweenLocations
{
public static double Calc(double Lat1,
double Long1, double Lat2, double Long2)
{
/*
The Haversine formula according to Dr. Math.
http://mathforum.org/library/drmath/view/51879.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2(sqrt(a), sqrt(1-a))
d = R * c
Where
* dlon is the change in longitude
* dlat is the change in latitude
* c is the great circle distance in Radians.
* R is the radius of a spherical Earth.
* The locations of the two points in
spherical coordinates (longitude and
latitude) are lon1,lat1 and lon2, lat2.
*/
double dDistance = Double.MinValue;
double dLat1InRad = Lat1 * (Math.PI / 180.0);
double dLong1InRad = Long1 * (Math.PI / 180.0);
double dLat2InRad = Lat2 * (Math.PI / 180.0);
double dLong2InRad = Long2 * (Math.PI / 180.0);
double dLongitude = dLong2InRad - dLong1InRad;
double dLatitude = dLat2InRad - dLat1InRad;
// Intermediate result a.
double a = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) +
Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) *
Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);
// Intermediate result c (great circle distance in Radians).
double c = 2.0 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1.0 - a));
// Distance.
// const Double kEarthRadiusMiles = 3956.0;
const Double kEarthRadiusKms = 6376.5;
dDistance = kEarthRadiusKms * c;
return dDistance;
}
public static double Calc(string NS1, double Lat1, double Lat1Min,
string EW1, double Long1, double Long1Min, string NS2,
double Lat2, double Lat2Min, string EW2,
double Long2, double Long2Min)
{
double NS1Sign = NS1.ToUpper() == "N" ? 1.0 : -1.0;
double EW1Sign = NS1.ToUpper() == "E" ? 1.0 : -1.0;
double NS2Sign = NS2.ToUpper() == "N" ? 1.0 : -1.0;
double EW2Sign = EW2.ToUpper() == "E" ? 1.0 : -1.0;
return (Calc(
(Lat1 + (Lat1Min / 60)) * NS1Sign,
(Long1 + (Long1Min / 60)) * EW1Sign,
(Lat2 + (Lat2Min / 60)) * NS2Sign,
(Long2 + (Long2Min / 60)) * EW2Sign
));
}
public static void Main(string[] args)
{
if (args.Length < 12)
{
System.Console.WriteLine("usage: DistanceBetweenLocations" +
" N 43 35.500 W 80 27.800 N 43 35.925 W 80 28.318");
return
}
System.Console.WriteLine(Calc(
args[0],
System.Double.Parse(args[1]),
System.Double.Parse(args[2]),
args[3],
System.Double.Parse(args[4]),
System.Double.Parse(args[5]),
args[6],
System.Double.Parse(args[7]),
System.Double.Parse(args[8]),
args[9],
System.Double.Parse(args[10]),
System.Double.Parse(args[11])));
}
}
어떠세요? 감이 오시나요? 다음은 위의 코드를 제가 C언어로 수정한 겁니다.
#include <math.h>
#define PI 3.14159265358979323846264338327950288419
double Calc(double Lat1, double Long1, double Lat2, double Long2)
{
/*
The Haversine formula according to Dr. Math.
http://mathforum.org/library/drmath/view/51879.html
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2(sqrt(a), sqrt(1-a))
d = R * c
Where
* dlon is the change in longitude
* dlat is the change in latitude
* c is the great circle distance in Radians.
* R is the radius of a spherical Earth.
* The locations of the two points in
spherical coordinates (longitude and
latitude) are lon1,lat1 and lon2, lat2.
*/
double dDistance;
const double kEarthRadiusKms = 6376.5;
double dLat1InRad = Lat1 * (PI / 180.0);
double dLong1InRad = Long1 * (PI / 180.0);
double dLat2InRad = Lat2 * (PI / 180.0);
double dLong2InRad = Long2 * (PI / 180.0);
double dLongitude = dLong2InRad - dLong1InRad;
double dLatitude = dLat2InRad - dLat1InRad;
// Intermediate result a.
double a = pow(sin(dLatitude / 2.0), 2.0) +
cos(dLat1InRad) * cos(dLat2InRad) *
pow(sin(dLongitude / 2.0), 2.0);
// Intermediate result c (great circle distance in Radians).
double c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a));
// Distance.
dDistance = kEarthRadiusKms * c;
return dDistance;
}
double Calc(char NS1, double Lat1, double Lat1Min,
char EW1, double Long1, double Long1Min, char NS2,
double Lat2, double Lat2Min, char EW2,
double Long2, double Long2Min)
{
double NS1Sign = (NS1 == 'N') ? 1.0 : -1.0;
double EW1Sign = (NS1 == 'E') ? 1.0 : -1.0;
double NS2Sign = (NS2 == 'N') ? 1.0 : -1.0;
double EW2Sign = (EW2 == 'E') ? 1.0 : -1.0;
return (Calc(
(Lat1 + (Lat1Min / 60)) * NS1Sign,
(Long1 + (Long1Min / 60)) * EW1Sign,
(Lat2 + (Lat2Min / 60)) * NS2Sign,
(Long2 + (Long2Min / 60)) * EW2Sign
));
}
int main()
{
//N 43 35.500 W 80 27.800 N 43 35.925 W 80 28.318
printf("Distance = %lf km \n",
Calc('N',43, 35.500, //Lat1
'W', 80, 27.800, //Lon1
'N', 43, 35.925, //Lat2
'W', 80, 28.318)); //Lon2
return 0;
}
자, 이게 전부입니다.
이 글은 제가 천문노트(http://astronote.org )에 이미 올린 글입니다.
참고글 : http://blog.jidolstar.com/217 (Flex로 만든 프로그램입니다.)
글쓴이 : 지돌스타(http://blog.jidolstar.com, http://astronote.org)