spheroid_direction()

double spheroid_direction	(	const GEOGRAPHIC_POINT *	r,
		const GEOGRAPHIC_POINT *	s,
		const SPHEROID *	spheroid
	)

Computes the direction of the geodesic joining two points on the spheroid.

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

r	- location of first point
s	- location of second point

Returns

azimuth of line joining r and s

Definition at line 287 of file lwspheroid.c.

{
        int i = 0;
        double lambda = s->lon - r->lon;
        double omf = 1 - spheroid->f;
        double u1 = atan(omf * tan(r->lat));
        double cos_u1 = cos(u1);
        double sin_u1 = sin(u1);
        double u2 = atan(omf * tan(s->lat));
        double cos_u2 = cos(u2);
        double sin_u2 = sin(u2);
 
        double omega = lambda;
        double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
        double sin_alpha, cos_alphasq, C, alphaFD;
        do
        {
                sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
                                POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
                sin_sigma = sqrt(sqr_sin_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 = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
                last_lambda = lambda;
                lambda = omega + (1.0 - C) * spheroid->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) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
 
        alphaFD = atan2((cos_u2 * sin(lambda)),
                        (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
        if (alphaFD < 0.0)
        {
                alphaFD = alphaFD + 2.0 * M_PI;
        }
        if (alphaFD > 2.0 * M_PI)
        {
                alphaFD = alphaFD - 2.0 * M_PI;
        }
        return alphaFD;
}

References SPHEROID::f, POW2, r, and s.

Referenced by lwgeom_azumith_spheroid(), and ptarray_area_spheroid().

Here is the caller graph for this function: