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;
    }
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)
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 .