cdm/src/main/java/ucar/unidata/geoloc/Bearing.java
public static Bearing calculateBearing(Earth e, double lat1, double lon1,
                                           double lat2, double lon2,
                                           Bearing result) {

        if (result == null) {
            result = new Bearing();
        }

        if ((lat1 == lat2) && (lon1 == lon2)) {
            result.distance    = 0;
            result.azimuth     = 0;
            result.backazimuth = 0;
            return result;
        }

        A = e.getMajor();
        F = e.getFlattening();
        R = 1.0 - F;

        // Algorithm from National Geodetic Survey, FORTRAN program "inverse,"
        // subroutine "INVER1," by L. PFEIFER and JOHN G. GERGEN.
        // http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
        // Conversion to JAVA from FORTRAN was made with as few changes as possible
        // to avoid errors made while recasting form, and to facilitate any future
        // comparisons between the original code and the altered version in Java.
        // Original documentation:
        // SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
        // MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
        // EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
        // STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
        // A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
        // F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERNECE ELLIPSOID
        // LATITUDES GLAT1 AND GLAT2
        // AND LONGITUDES GLON1 AND GLON2 ARE IN RADIANS POSITIVE NORTH AND EAST
        // FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH
        //
        // Reference ellipsoid is the WGS-84 ellipsoid.
        // See http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
        // FAZ is forward azimuth in radians from pt1 to pt2;
        // BAZ is backward azimuth from point 2 to 1;
        // S is distance in meters.
        //
        // Conversion to JAVA from FORTRAN was made with as few changes as possible
        // to avoid errors made while recasting form, and to facilitate any future
        // comparisons between the original code and the altered version in Java.
        //
        //IMPLICIT REAL*8 (A-H,O-Z)
        //  COMMON/CONST/PI,RAD
        //  COMMON/ELIPSOID/A,F
        double GLAT1 = rad * lat1;
        double GLAT2 = rad * lat2;
        double TU1   = R * Math.sin(GLAT1) / Math.cos(GLAT1);
        double TU2   = R * Math.sin(GLAT2) / Math.cos(GLAT2);
        double CU1   = 1. / Math.sqrt(TU1 * TU1 + 1.);
        double SU1   = CU1 * TU1;
        double CU2   = 1. / Math.sqrt(TU2 * TU2 + 1.);
        double S     = CU1 * CU2;
        double BAZ   = S * TU2;
        double FAZ   = BAZ * TU1;
        double GLON1 = rad * lon1;
        double GLON2 = rad * lon2;
        double X     = GLON2 - GLON1;
        double D, SX, CX, SY, CY, Y, SA, C2A, CZ, E, C;
        int    loopCnt = 0;
        do {
            loopCnt++;
            //Check for an infinite loop
            if (loopCnt > 1000) {
                throw new IllegalArgumentException(
                    "Too many iterations calculating bearing:" + lat1 + " "
                    + lon1 + " " + lat2 + " " + lon2);
            }
            SX  = Math.sin(X);
            CX  = Math.cos(X);
            TU1 = CU2 * SX;
            TU2 = BAZ - SU1 * CU2 * CX;
            SY  = Math.sqrt(TU1 * TU1 + TU2 * TU2);
            CY  = S * CX + FAZ;
            Y   = Math.atan2(SY, CY);
            SA  = S * SX / SY;
            C2A = -SA * SA + 1.;
            CZ  = FAZ + FAZ;
            if (C2A > 0.) {
                CZ = -CZ / C2A + CY;
            }
            E = CZ * CZ * 2. - 1.;
            C = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
            D = X;
            X = ((E * CY * C + CZ) * SY * C + Y) * SA;
            X = (1. - C) * X * F + GLON2 - GLON1;
            //IF(DABS(D-X).GT.EPS) GO TO 100
        } while (Math.abs(D - X) > EPS);

        if (loopCnt > maxLoopCnt) {
            maxLoopCnt = loopCnt;
            //        System.err.println("loopCnt:" + loopCnt);
        }

        FAZ = Math.atan2(TU1, TU2);
        BAZ = Math.atan2(CU1 * SX, BAZ * CX - SU1 * CU2) + Math.PI;
        X   = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
        X   = (X - 2.) / X;
        C   = 1. - X;
        C   = (X * X / 4. + 1.) / C;
        D   = (0.375 * X * X - 1.) * X;
        X   = E * CY;
        S   = 1. - E - E;
        S = ((((SY * SY * 4. - 3.) * S * CZ * D / 6. - X) * D / 4. + CZ) * SY
             * D + Y) * C * A * R;

        result.distance = S / 1000.0;  // meters to km
        result.azimuth  = FAZ * deg;   // radians to degrees

        if (result.azimuth < 0.0) {
            result.azimuth += 360.0;  // reset azs from -180 to 180 to 0 to 360
        }

        result.backazimuth = BAZ * deg;  // radians to degrees; already in 0 to 360 range

        return result;
    }
public static Bearing calculateBearing(Earth e, double lat1, double lon1,
                                           double lat2, double lon2,
                                           Bearing result) {

        if (result == null) {
            result = new Bearing();
        }

        if ((lat1 == lat2) && (lon1 == lon2)) {
            result.distance    = 0;
            result.azimuth     = 0;
            result.backazimuth = 0;
            return result;
        }

        A = e.getMajor();
        F = e.getFlattening();
        R = 1.0 - F;

        // Algorithm from National Geodetic Survey, FORTRAN program "inverse,"
        // subroutine "INVER1," by L. PFEIFER and JOHN G. GERGEN.
        // http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
        // Conversion to JAVA from FORTRAN was made with as few changes as possible
        // to avoid errors made while recasting form, and to facilitate any future
        // comparisons between the original code and the altered version in Java.
        // Original documentation:
        // SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
        // MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
        // EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
        // STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
        // A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
        // F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERNECE ELLIPSOID
        // LATITUDES GLAT1 AND GLAT2
        // AND LONGITUDES GLON1 AND GLON2 ARE IN RADIANS POSITIVE NORTH AND EAST
        // FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH
        //
        // Reference ellipsoid is the WGS-84 ellipsoid.
        // See http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
        // FAZ is forward azimuth in radians from pt1 to pt2;
        // BAZ is backward azimuth from point 2 to 1;
        // S is distance in meters.
        //
        // Conversion to JAVA from FORTRAN was made with as few changes as possible
        // to avoid errors made while recasting form, and to facilitate any future
        // comparisons between the original code and the altered version in Java.
        //
        //IMPLICIT REAL*8 (A-H,O-Z)
        //  COMMON/CONST/PI,RAD
        //  COMMON/ELIPSOID/A,F
        double GLAT1 = rad * lat1;
        double GLAT2 = rad * lat2;
        double TU1   = R * Math.sin(GLAT1) / Math.cos(GLAT1);
        double TU2   = R * Math.sin(GLAT2) / Math.cos(GLAT2);
        double CU1   = 1. / Math.sqrt(TU1 * TU1 + 1.);
        double SU1   = CU1 * TU1;
        double CU2   = 1. / Math.sqrt(TU2 * TU2 + 1.);
        double S     = CU1 * CU2;
        double BAZ   = S * TU2;
        double FAZ   = BAZ * TU1;
        double GLON1 = rad * lon1;
        double GLON2 = rad * lon2;
        double X     = GLON2 - GLON1;
        double D, SX, CX, SY, CY, Y, SA, C2A, CZ, E, C;
        int    loopCnt = 0;
        do {
            loopCnt++;
            //Check for an infinite loop
            if (loopCnt > 1000) {
                throw new IllegalArgumentException(
                    "Too many iterations calculating bearing:" + lat1 + " "
                    + lon1 + " " + lat2 + " " + lon2);
            }
            SX  = Math.sin(X);
            CX  = Math.cos(X);
            TU1 = CU2 * SX;
            TU2 = BAZ - SU1 * CU2 * CX;
            SY  = Math.sqrt(TU1 * TU1 + TU2 * TU2);
            CY  = S * CX + FAZ;
            Y   = Math.atan2(SY, CY);
            SA  = S * SX / SY;
            C2A = -SA * SA + 1.;
            CZ  = FAZ + FAZ;
            if (C2A > 0.) {
                CZ = -CZ / C2A + CY;
            }
            E = CZ * CZ * 2. - 1.;
            C = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
            D = X;
            X = ((E * CY * C + CZ) * SY * C + Y) * SA;
            X = (1. - C) * X * F + GLON2 - GLON1;
            //IF(DABS(D-X).GT.EPS) GO TO 100
        } while (Math.abs(D - X) > EPS);

        if (loopCnt > maxLoopCnt) {
            maxLoopCnt = loopCnt;
            //        System.err.println("loopCnt:" + loopCnt);
        }

        FAZ = Math.atan2(TU1, TU2);
        BAZ = Math.atan2(CU1 * SX, BAZ * CX - SU1 * CU2) + Math.PI;
        X   = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
        X   = (X - 2.) / X;
        C   = 1. - X;
        C   = (X * X / 4. + 1.) / C;
        D   = (0.375 * X * X - 1.) * X;
        X   = E * CY;
        S   = 1. - E - E;
        S = ((((SY * SY * 4. - 3.) * S * CZ * D / 6. - X) * D / 4. + CZ) * SY
             * D + Y) * C * A * R;

        result.distance = S / 1000.0;  // meters to km
        result.azimuth  = FAZ * deg;   // radians to degrees

        if (result.azimuth < 0.0) {
            result.azimuth += 360.0;  // reset azs from -180 to 180 to 0 to 360
        }

        result.backazimuth = BAZ * deg;  // radians to degrees; already in 0 to 360 range

        return result;
    }
Computes distance (in km), azimuth (degrees clockwise positive
from North, 0 to 360), and back azimuth (degrees clockwise positive
from North, 0 to 360), from latitude-longituide point pt1 to
latitude-longituide pt2.<p>
Algorithm from U.S. National Geodetic Survey, FORTRAN program "inverse,"
subroutine "INVER1," by L. PFEIFER and JOHN G. GERGEN.
See http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
<P>Original documentation:
<br>SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
<br>MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
<br>EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
<br>STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
</P>
Reference ellipsoid is the WGS-84 ellipsoid.
<br>See http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
<p/>
Requires close to 1.4 E-5 seconds wall clock time per call
on a 550 MHz Pentium with Linux 7.2.

@param e      Earth object (defines radius and flattening)
@param lat1   Lat of point 1
@param lon1   Lon of point 1
@param lat2   Lat of point 2
@param lon2   Lon of point 2
@param result put result here, or null to allocate
@return a Bearing object with distance (in km), azimuth from
pt1 to pt2 (degrees, 0 = north, clockwise positive)