spheroid_project()

int spheroid_project	(	const GEOGRAPHIC_POINT *	r,
		const SPHEROID *	spheroid,
		double	distance,
		double	azimuth,
		GEOGRAPHIC_POINT *	g
	)

Given a location, an azimuth and a distance, computes the location of the projected point.

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

Parameters

r	- location of first point.
distance	- distance in meters.
azimuth	- azimuth in radians.

Returns

s - location of projected point.

Definition at line 364 of file lwspheroid.c.

{
        double omf = 1 - spheroid->f;
        double tan_u1 = omf * tan(r->lat);
        double u1 = atan(tan_u1);
        double sigma, last_sigma, delta_sigma, two_sigma_m;
        double sigma1, sin_alpha, alpha, cos_alphasq;
        double u2, A, B;
        double lat2, lambda, lambda2, C, omega;
        int i = 0;
 
        if (azimuth < 0.0)
        {
                azimuth = azimuth + M_PI * 2.0;
        }
        if (azimuth > (M_PI * 2.0))
        {
                azimuth = azimuth - M_PI * 2.0;
        }
 
        sigma1 = atan2(tan_u1, cos(azimuth));
        sin_alpha = cos(u1) * sin(azimuth);
        alpha = asin(sin_alpha);
        cos_alphasq = 1.0 - POW2(sin_alpha);
 
        u2 = spheroid_mu2(alpha, spheroid);
        A = spheroid_big_a(u2);
        B = spheroid_big_b(u2);
 
        sigma = (distance / (spheroid->b * A));
        do
        {
                two_sigma_m = 2.0 * sigma1 + sigma;
                delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
                last_sigma = sigma;
                sigma = (distance / (spheroid->b * A)) + delta_sigma;
                i++;
        }
        while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
 
        lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
                      cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
                                                 POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
                                                      cos(azimuth)))));
        lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
                       sin(u1) * sin(sigma) * cos(azimuth)));
        C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
        omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
                (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
        lambda2 = r->lon + omega;
        g->lat = lat2;
        g->lon = lambda2;
        return LW_SUCCESS;
}

References SPHEROID::b, distance(), SPHEROID::f, GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_SUCCESS, POW2, r, spheroid_big_a(), spheroid_big_b(), and spheroid_mu2().

Referenced by lwgeom_project_spheroid(), and ptarray_area_spheroid().

Here is the call graph for this function:

Here is the caller graph for this function: