- double *course, double *dist ) {
- // d = 2*asin(sqrt((sin((lat1-lat2)/2))^2 +
- // cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))
- double tmp1 = sin( (start.y() - dest.y()) / 2 );
- double tmp2 = sin( (start.x() - dest.x()) / 2 );
- double d = 2.0 * asin( sqrt( tmp1 * tmp1 +
- cos(start.y()) * cos(dest.y()) * tmp2 * tmp2));
-
- // We obtain the initial course, tc1, (at point 1) from point 1 to
- // point 2 by the following. The formula fails if the initial
- // point is a pole. We can special case this with:
- //
- // IF (cos(lat1) < EPS) // EPS a small number ~ machine precision
- // IF (lat1 > 0)
- // tc1= pi // starting from N pole
- // ELSE
- // tc1= 0 // starting from S pole
- // ENDIF
- // ENDIF
- //
- // For starting points other than the poles:
- //
- // IF sin(lon2-lon1)<0
- // tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
- // ELSE
- // tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
- // ENDIF
-
- double tc1;
-
- if ( cos(start.y()) < SG_EPSILON ) {
- // EPS a small number ~ machine precision
- if ( start.y() > 0 ) {
- tc1 = SGD_PI; // starting from N pole
- } else {
- tc1 = 0; // starting from S pole
- }
- }
-
- // For starting points other than the poles:
-
- double tmp3 = sin(d)*cos(start.y());
- double tmp4 = sin(dest.y())-sin(start.y())*cos(d);
- double tmp5 = acos(tmp4/tmp3);
- if ( sin( dest.x() - start.x() ) < 0 ) {
- tc1 = tmp5;
- } else {
- tc1 = 2 * SGD_PI - tmp5;
- }
-
- *course = tc1;
- *dist = d * SG_RAD_TO_NM * SG_NM_TO_METER;
+ double *course, double *dist )
+{
+ SGGeoc gs = start.toSGGeoc();
+ SGGeoc gd = dest.toSGGeoc();
+ *course = SGGeoc::courseRad(gs, gd);
+ *dist = SGGeoc::distanceM(gs, gd);