spheroid_distance()

double spheroid_distance	(	const GEOGRAPHIC_POINT *	a,
		const GEOGRAPHIC_POINT *	b,
		const SPHEROID *	spheroid
	)

Computes the shortest distance along the surface of the spheroid between two points.

Based on Vincenty's formula for the geodetic inverse problem as described in "Geocentric Datum of Australia Technical Manual", Chapter 4. Tested against: http://mascot.gdbc.gov.bc.ca/mascot/util1a.html and http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp

Parameters

a	- location of first point.
b	- location of second point.
s	- spheroid to calculate on

Returns

spheroidal distance between a and b in spheroid units.

Definition at line 186 of file lwspheroid.c.

References SPHEROID::a, SPHEROID::b, distance(), SPHEROID::f, geographic_point_equals(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, lwerror(), POW2, SPHEROID::radius, sphere_distance(), spheroid_big_a(), spheroid_big_b(), and spheroid_mu2().

Referenced by circ_tree_distance_tree(), ptarray_area_spheroid(), ptarray_distance_spheroid(), ptarray_length_spheroid(), spheroid_big_b(), and test_spheroid_distance().

{
        double lambda = (b->lon - a->lon);
        double f = spheroid->f;
        double omf = 1 - spheroid->f;
        double u1, u2;
        double cos_u1, cos_u2;
        double sin_u1, sin_u2;
        double big_a, big_b, delta_sigma;
        double alpha, sin_alpha, cos_alphasq, c;
        double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
        double cos_lambda, sin_lambda;
        double distance;
        int i = 0;

        /* Same point => zero distance */
        if ( geographic_point_equals(a, b) )
        {
                return 0.0;
        }

        u1 = atan(omf * tan(a->lat));
        cos_u1 = cos(u1);
        sin_u1 = sin(u1);
        u2 = atan(omf * tan(b->lat));
        cos_u2 = cos(u2);
        sin_u2 = sin(u2);

        omega = lambda;
        do
        {
                cos_lambda = cos(lambda);
                sin_lambda = sin(lambda);
                sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
                               POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
                sin_sigma = sqrt(sqrsin_sigma);
                cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
                sigma = atan2(sin_sigma, cos_sigma);
                sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);

                /* Numerical stability issue, ensure asin is not NaN */
                if ( sin_alpha > 1.0 )
                        alpha = M_PI_2;
                else if ( sin_alpha < -1.0 )
                        alpha = -1.0 * M_PI_2;
                else
                        alpha = asin(sin_alpha);

                cos_alphasq = POW2(cos(alpha));
                cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);

                /* Numerical stability issue, cos2 is in range */
                if ( cos2_sigma_m > 1.0 )
                        cos2_sigma_m = 1.0;
                if ( cos2_sigma_m < -1.0 )
                        cos2_sigma_m = -1.0;

                c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
                last_lambda = lambda;
                lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
                         (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
                i++;
        }
        while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );

        u2 = spheroid_mu2(alpha, spheroid);
        big_a = spheroid_big_a(u2);
        big_b = spheroid_big_b(u2);
        delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
                                           (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));

        distance = spheroid->b * big_a * (sigma - delta_sigma);

        /* Algorithm failure, distance == NaN, fallback to sphere */
        if ( distance != distance )
        {
                lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a->lat, a->lon, b->lat, b->lon, spheroid->a, spheroid->b);
                return spheroid->radius * sphere_distance(a, b);
        }

        return distance;
}

Here is the call graph for this function:

Here is the caller graph for this function: