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