cdm/src/main/java/ucar/unidata/geoloc/Bearing.java
public static LatLonPointImpl findPoint(Earth e, double lat1,
                                            double lon1, double az,
                                            double dist,
                                            LatLonPointImpl result) {
        if (result == null) {
            result = new LatLonPointImpl();
        }

        if ((dist == 0)) {
            result.setLatitude(lat1);
            result.setLongitude(lon1);
            return result;
        }

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

        // Algorithm from National Geodetic Survey, FORTRAN program "forward,"
        // subroutine "DIRCT1," by stephen j. frakes.
        // 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:
        //   SUBROUTINE DIRCT1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S)
        //
        //   SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY
        //   MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
        //   EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
        //
        //   A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
        //   F IS THE FLATTENING OF THE REFERENCE ELLIPSOID
        //   LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
        //   AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH
        //   GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A
        //
        //   PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75
        //   MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608
        //

        if (az < 0.0) {
            az += 360.0;  // reset azs from -180 to 180 to 0 to 360
        }
        double FAZ   = az * rad;
        double GLAT1 = lat1 * rad;
        double GLON1 = lon1 * rad;
        double S     = dist * 1000.;  // convert to meters
        double TU    = R * Math.sin(GLAT1) / Math.cos(GLAT1);
        double SF    = Math.sin(FAZ);
        double CF    = Math.cos(FAZ);
        double BAZ   = 0.;
        if (CF != 0) {
            BAZ = Math.atan2(TU, CF) * 2;
        }
        double CU  = 1. / Math.sqrt(TU * TU + 1.);
        double SU  = TU * CU;
        double SA  = CU * SF;
        double C2A = -SA * SA + 1.;
        double X   = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
        X = (X - 2.) / X;
        double C = 1. - X;
        C = (X * X / 4. + 1) / C;
        double D = (0.375 * X * X - 1.) * X;
        TU = S / R / A / C;
        double Y = TU;
        double SY, CY, CZ, E, GLAT2, GLON2;
        do {
            SY = Math.sin(Y);
            CY = Math.cos(Y);
            CZ = Math.cos(BAZ + Y);
            E  = CZ * CZ * 2. - 1.;
            C  = Y;
            X  = E * CY;
            Y  = E + E - 1.;
            Y = (((SY * SY * 4. - 3.) * Y * CZ * D / 6. + X) * D / 4. - CZ)
                * SY * D + TU;
        } while (Math.abs(Y - C) > EPS);
        BAZ   = CU * CY * CF - SU * SY;
        C     = R * Math.sqrt(SA * SA + BAZ * BAZ);
        D     = SU * CY + CU * SY * CF;
        GLAT2 = Math.atan2(D, C);
        C     = CU * CY - SU * SY * CF;
        X     = Math.atan2(SY * SF, C);
        C     = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
        D     = ((E * CY * C + CZ) * SY * C + Y) * SA;
        GLON2 = GLON1 + X - (1. - C) * D * F;
        BAZ   = (Math.atan2(SA, BAZ) + Math.PI) * deg;
        result.setLatitude(GLAT2 * deg);
        result.setLongitude(GLON2 * deg);
        return result;
    }
public static LatLonPointImpl findPoint(Earth e, double lat1,
                                            double lon1, double az,
                                            double dist,
                                            LatLonPointImpl result) {
        if (result == null) {
            result = new LatLonPointImpl();
        }

        if ((dist == 0)) {
            result.setLatitude(lat1);
            result.setLongitude(lon1);
            return result;
        }

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

        // Algorithm from National Geodetic Survey, FORTRAN program "forward,"
        // subroutine "DIRCT1," by stephen j. frakes.
        // 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:
        //   SUBROUTINE DIRCT1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S)
        //
        //   SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY
        //   MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
        //   EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
        //
        //   A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
        //   F IS THE FLATTENING OF THE REFERENCE ELLIPSOID
        //   LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
        //   AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH
        //   GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A
        //
        //   PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75
        //   MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608
        //

        if (az < 0.0) {
            az += 360.0;  // reset azs from -180 to 180 to 0 to 360
        }
        double FAZ   = az * rad;
        double GLAT1 = lat1 * rad;
        double GLON1 = lon1 * rad;
        double S     = dist * 1000.;  // convert to meters
        double TU    = R * Math.sin(GLAT1) / Math.cos(GLAT1);
        double SF    = Math.sin(FAZ);
        double CF    = Math.cos(FAZ);
        double BAZ   = 0.;
        if (CF != 0) {
            BAZ = Math.atan2(TU, CF) * 2;
        }
        double CU  = 1. / Math.sqrt(TU * TU + 1.);
        double SU  = TU * CU;
        double SA  = CU * SF;
        double C2A = -SA * SA + 1.;
        double X   = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
        X = (X - 2.) / X;
        double C = 1. - X;
        C = (X * X / 4. + 1) / C;
        double D = (0.375 * X * X - 1.) * X;
        TU = S / R / A / C;
        double Y = TU;
        double SY, CY, CZ, E, GLAT2, GLON2;
        do {
            SY = Math.sin(Y);
            CY = Math.cos(Y);
            CZ = Math.cos(BAZ + Y);
            E  = CZ * CZ * 2. - 1.;
            C  = Y;
            X  = E * CY;
            Y  = E + E - 1.;
            Y = (((SY * SY * 4. - 3.) * Y * CZ * D / 6. + X) * D / 4. - CZ)
                * SY * D + TU;
        } while (Math.abs(Y - C) > EPS);
        BAZ   = CU * CY * CF - SU * SY;
        C     = R * Math.sqrt(SA * SA + BAZ * BAZ);
        D     = SU * CY + CU * SY * CF;
        GLAT2 = Math.atan2(D, C);
        C     = CU * CY - SU * SY * CF;
        X     = Math.atan2(SY * SF, C);
        C     = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
        D     = ((E * CY * C + CZ) * SY * C + Y) * SA;
        GLON2 = GLON1 + X - (1. - C) * D * F;
        BAZ   = (Math.atan2(SA, BAZ) + Math.PI) * deg;
        result.setLatitude(GLAT2 * deg);
        result.setLongitude(GLON2 * deg);
        return result;
    }
public static LatLonPointImpl findPoint ( Earth e , double lat1 , double lon1 , double az , double dist , LatLonPointImpl result ) { if ( result == null ) { result = new LatLonPointImpl ( ) ; } if ( ( dist == 0 ) ) { result . setLatitude ( lat1 ) ; result . setLongitude ( lon1 ) ; return result ; } A = e . getMajor ( ) ; F = e . getFlattening ( ) ; R = 1.0 - F ; // Algorithm from National Geodetic Survey, FORTRAN program "forward,"
 // subroutine "DIRCT1," by stephen j. frakes.
 // 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:
 //   SUBROUTINE DIRCT1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S)
 //
 //   SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY
 //   MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
 //   EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
 //
 //   A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
 //   F IS THE FLATTENING OF THE REFERENCE ELLIPSOID
 //   LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
 //   AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH
 //   GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A
 //
 //   PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75
 //   MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608
 //
 if ( az < 0.0 ) { az += 360.0 ; // reset azs from -180 to 180 to 0 to 360
 } double FAZ = az * rad ; double GLAT1 = lat1 * rad ; double GLON1 = lon1 * rad ; double S = dist * 1000. ; // convert to meters
 double TU = R * Math . sin ( GLAT1 ) / Math . cos ( GLAT1 ) ; double SF = Math . sin ( FAZ ) ; double CF = Math . cos ( FAZ ) ; double BAZ = 0. ; if ( CF != 0 ) { BAZ = Math . atan2 ( TU , CF ) * 2 ; } double CU = 1. / Math . sqrt ( TU * TU + 1. ) ; double SU = TU * CU ; double SA = CU * SF ; double C2A = - SA * SA + 1. ; double X = Math . sqrt ( ( 1. / R / R - 1. ) * C2A + 1. ) + 1. ; X = ( X - 2. ) / X ; double C = 1. - X ; C = ( X * X / 4. + 1 ) / C ; double D = ( 0.375 * X * X - 1. ) * X ; TU = S / R / A / C ; double Y = TU ; double SY , CY , CZ , E , GLAT2 , GLON2 ; do { SY = Math . sin ( Y ) ; CY = Math . cos ( Y ) ; CZ = Math . cos ( BAZ + Y ) ; E = CZ * CZ * 2. - 1. ; C = Y ; X = E * CY ; Y = E + E - 1. ; Y = ( ( ( SY * SY * 4. - 3. ) * Y * CZ * D / 6. + X ) * D / 4. - CZ ) * SY * D + TU ; } while ( Math . abs ( Y - C ) > EPS ) ; BAZ = CU * CY * CF - SU * SY ; C = R * Math . sqrt ( SA * SA + BAZ * BAZ ) ; D = SU * CY + CU * SY * CF ; GLAT2 = Math . atan2 ( D , C ) ; C = CU * CY - SU * SY * CF ; X = Math . atan2 ( SY * SF , C ) ; C = ( ( - 3. * C2A + 4. ) * F + 4. ) * C2A * F / 16. ; D = ( ( E * CY * C + CZ ) * SY * C + Y ) * SA ; GLON2 = GLON1 + X - ( 1. - C ) * D * F ; BAZ = ( Math . atan2 ( SA , BAZ ) + Math . PI ) * deg ; result . setLatitude ( GLAT2 * deg ) ; result . setLongitude ( GLON2 * deg ) ; return result ; }
Calculate a position given an azimuth and distance from
another point.
<p/>
<p/>
Algorithm from National Geodetic Survey, FORTRAN program "forward,"
subroutine "DIRCT1," by stephen j. frakes.
http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
<p>Original documentation:
<pre>
SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY
MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
</pre>

@param e      Earth object (defines radius and flattening)
@param lat1   latitude of starting point
@param lon1   longitude of starting point
@param az     forward azimuth (degrees)
@param dist   distance from the point (km)
@param result Object to use if non-null
@return the position as a LatLonPointImpl
Calculate a position given an azimuth and distance from another point . <p / > <p / > Algorithm from National Geodetic Survey FORTRAN program forward subroutine DIRCT1 by stephen j . frakes . http : // www . ngs . noaa . gov / TOOLS / Inv_Fwd / Inv_Fwd . html <p > Original documentation : <pre > SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T . VINCENTY MODIFIED RAINSFORD S METHOD WITH HELMERT S ELLIPTICAL TERMS EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL < / pre >