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 282 of file lwspheroid.c.

References SPHEROID::f, GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, and POW2.

Referenced by lwgeom_azumith_spheroid(), ptarray_area_spheroid(), and spheroid_big_b().

{
        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;
}

Here is the caller graph for this function: