distance_ellipse_calculation()

double distance_ellipse_calculation	(	double	lat1,
		double	long1,
		double	lat2,
		double	long2,
		SPHEROID *	sphere
	)

Definition at line 246 of file lwgeom_spheroid.c.

References SPHEROID::b, bigA(), bigB(), deltaLongitude(), SPHEROID::f, LWGEOM_length2d_ellipsoid(), mu2(), and PG_FUNCTION_INFO_V1().

Referenced by distance_ellipse().

{
        /*
         * Code is taken from David Skea
         * Geographic Data BC, Province of British Columbia, Canada.
         *  Thanks to GDBC and David Skea for allowing this to be
         *  put in PostGIS.
         */

        double L1,L2,sinU1,sinU2,cosU1,cosU2;
        double dl,dl1,dl2,dl3,cosdl1,sindl1;
        double cosSigma,sigma,azimuthEQ,tsm;
        double u2,A,B;
        double dsigma;

        double TEMP;

        int iterations;


        L1 = atan((1.0 - sphere->f ) * tan( lat1) );
        L2 = atan((1.0 - sphere->f ) * tan( lat2) );
        sinU1 = sin(L1);
        sinU2 = sin(L2);
        cosU1 = cos(L1);
        cosU2 = cos(L2);

        dl = long2- long1;
        dl1 = dl;
        cosdl1 = cos(dl);
        sindl1 = sin(dl);
        iterations = 0;
        do
        {
                cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
                sigma = acos(cosSigma);
                azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma));

                /*
                 * Patch from Patrica Tozer to handle minor
                 * mathematical stability problem
                 */
                TEMP = cosSigma - (2.0 * sinU1 * sinU2)/
                       (cos(azimuthEQ)*cos(azimuthEQ));
                if (TEMP > 1)
                {
                        TEMP = 1;
                }
                else if (TEMP < -1)
                {
                        TEMP = -1;
                }
                tsm = acos(TEMP);


                /* (old code?)
                tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
                */

                dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
                dl3 = dl1 - (dl + dl2);
                dl1 = dl + dl2;
                cosdl1 = cos(dl1);
                sindl1 = sin(dl1);
                iterations++;
        }
        while ( (iterations<999) && (fabs(dl3) > 1.0e-032));

        /* compute expansions A and B */
        u2 = mu2(azimuthEQ,sphere);
        A = bigA(u2);
        B = bigB(u2);

        /* compute length of geodesic */
        dsigma = B * sin(sigma) * (cos(tsm) +
                                   (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0);
        return sphere->b * (A * (sigma - dsigma));
}

Here is the call graph for this function:

Here is the caller graph for this function: