/*
 * $Id: Bearing.java,v 1.19 2004/12/10 15:07:50 dmurray Exp $
 *
 * Copyright © 1997-2004 Unidata Program Center/University Corporation for
 * Atmospheric Research, P.O. Box 3000, Boulder, CO 80307,
 * support@unidata.ucar.edu.
 *
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at
 * your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
 * General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

package ucar.unidata.geoloc;



import java.lang.Math;


/**
 * Computes the distance, azimuth, and back azimuth between
 * two lat-lon positions on the Earth's surface. Reference ellipsoid is the WGS-84.
 *
 * @author Unidata Development Team
 * @version $Id: Bearing.java,v 1.19 2004/12/10 15:07:50 dmurray Exp $
 */
public class Bearing {

    /** the azimuth, degrees, 0 = north, clockwise positive */
    private double azimuth;

    /** the back azimuth, degrees, 0 = north, clockwise positive */
    private double backazimuth;

    /** separation in kilometers */
    private double distance;

    /**
     * Get the azimuth in degrees, 0 = north, clockwise positive
     *
     * @return azimuth in degrees
     */
    public double getAngle() {
        return azimuth;
    }

    /**
     * Get the back azimuth in degrees, 0 = north, clockwise positive
     * @return back azimuth in degrees
     */
    public double getBackAzimuth() {
        return backazimuth;
    }

    /**
     * Get the distance in kilometers
     * @return distance in km
     */
    public double getDistance() {
        return distance;
    }

    /**
     * Calculate the bearing between the 2 points.
     * See calculateBearing below.
     *
     * @param pt1 Point 1
     * @param pt2 Point 2
     * @param result Object to use if non-null
     *
     * @return The bearing
     */
    public static Bearing calculateBearing(LatLonPoint pt1, LatLonPoint pt2,
                                           Bearing result) {
        return calculateBearing(pt1.getLatitude(), pt1.getLongitude(),
                                pt2.getLatitude(), pt2.getLongitude(),
                                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
     *
     * Requires close to 1.4 E-5 seconds wall clock time per call
     * on a 550 MHz Pentium with Linux 7.2.
     *
     * @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)
     */
    public static Bearing calculateBearing(double lat1, double lon1,
                                           double lat2, double lon2,
                                           Bearing result) {

        // to convert degrees to radians, multiply by:
        final double rad = Math.toRadians(1.0);
        // to convert radians to degrees:
        final double deg = Math.toDegrees(1.0);

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

        // 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 A = 6378137.0;  // in meters
        double F = 1.0 / 298.257223563;
        // DATA EPS/0.5D-13/
        double EPS   = 0.5E-13;
        double R     = 1.0 - 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;
        do {
            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);

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


    /*
     * This method is for same use as calculateBearing, but has much simpler calculations
     * by assuming a spherical earth. It is actually slower than
     * "calculateBearing" code, probably due to having more trig function calls.
     * It is less accurate, too.
     * Errors are on the order of 1/300 or less. This code
     * saved here only as a warning to future programmers thinking of this approach.
     *
     * Requires close to 2.0 E-5 seconds wall clock time per call
     * on a 550 MHz Pentium with Linux 7.2.
     *
     * public static Bearing calculateBearingAlternate
     *   (LatLonPoint pt1, LatLonPoint pt2, Bearing result) {
     *
     * // to convert degrees to radians, multiply by:
     * final double rad = Math.toRadians(1.0);
     * // to convert radians to degrees:
     * final double deg = Math.toDegrees(1.0);
     *
     * if (result == null)
     *   result = new Bearing();
     *
     * double R = 6371008.7;  // mean earth radius in meters; WGS 84 definition
     * double GLAT1 = rad*(pt1.getLatitude());
     * double GLAT2 = rad*(pt2.getLatitude());
     * double GLON1 = rad*(pt1.getLongitude());
     * double GLON2 = rad*(pt2.getLongitude());
     *
     * // great circle angular separation in radians
     * double alpha = Math.acos( Math.sin(GLAT1)*Math.sin(GLAT2)
     *                          +Math.cos(GLAT1)*Math.cos(GLAT2)*Math.cos(GLON1-GLON2) );
     * // great circle distance in meters
     * double gcd = R * alpha;
     *
     * result.distance = gcd / 1000.0;      // meters to km
     *
     * // forward azimuth from point 1 to 2 in radians
     * double s2 = rad*(90.0-pt2.getLatitude());
     * double FAZ = Math.asin(Math.sin(s2)*Math.sin(GLON2-GLON1) / Math.sin(alpha));
     *
     * result.azimuth = FAZ * deg;        // radians to degrees
     * if (result.azimuth < 0.0)
     * result.azimuth += 360.0;       // reset az from -180 to 180 to 0 to 360
     *
     * // back azimuth from point 2 to 1 in radians
     * double s1 = rad*(90.0-pt1.getLatitude());
     * double BAZ = Math.asin(Math.sin(s1)*Math.sin(GLON1-GLON2) / Math.sin(alpha));
     *
     * result.backazimuth = BAZ * deg;
     * if (result.backazimuth < 0.0)
     *     result.backazimuth += 360.0;   // reset backaz from -180 to 180 to 0 to 360
     *
     * return result;
     * }
     */


    /**
     * Nice format.
     *
     * @return  return a nice format of this Bearing
     */
    public String toString() {
        StringBuffer buf = new StringBuffer();
        buf.append("Azimuth: ");
        buf.append(azimuth);
        buf.append("Back azimuth: ");
        buf.append(backazimuth);
        buf.append("Distance: ");
        buf.append(distance);
        return buf.toString();
    }

}

/* Change History:
   $Log: Bearing.java,v $
   Revision 1.19  2004/12/10 15:07:50  dmurray
   Jindent John's changes

   Revision 1.18  2004/09/22 21:21:55  caron
   merge; remove visad dependednt method


   Revision 1.15  2004/07/30 17:22:19  dmurray
   Jindent and doclint

   Revision 1.14  2004/07/30 16:24:39  dmurray
   Jindent and javadoc

   Revision 1.13  2004/02/27 21:21:27  jeffmc
   Lots of javadoc warning fixes

   Revision 1.12  2004/01/29 17:34:57  jeffmc
   A big sweeping checkin after a big sweeping reformatting
   using the new jindent.

   jindent adds in javadoc templates and reformats existing javadocs. In the new javadoc
   templates there is a '_more_' to remind us to fill these in.

   Revision 1.11  2003/04/08 15:59:05  caron
   rework for nc2 framework

   Revision 1.10  2002/11/13 16:40:47  wier
   revise spherical earth code

   Revision 1.9  2002/11/12 23:53:30  wier
   new method calculateBearingFast to compute distance and azimuth between two lat,long points faster than more accurate method

   Revision 1.8  2002/09/06 15:37:25  wier
   add back azimuth to calculation results

   Revision 1.7  2002/09/05 22:11:58  wier
   improve calculations to use ellipsoidal earth

   Revision 1.6  2002/05/10 13:59:10  dmurray
   added a toString method

   Revision 1.5  2000/08/18 04:15:15  russ
   Licensed under GNU LGPL.

   Revision 1.4  1999/12/16 22:57:16  caron
   gridded data viewer checkin

   Revision 1.3  1999/06/03 01:43:48  caron
   remove the damn controlMs

   Revision 1.2  1999/06/03 01:26:13  caron
   another reorg

   Revision 1.1.1.1  1999/05/21 17:33:51  caron
   startAgain

# Revision 1.5  1999/03/16  16:56:46  caron
# fix StationModel editing; add TopLevel
#
# Revision 1.4  1999/03/03  19:58:07  caron
# more java2D changes
#
# Revision 1.3  1999/02/15  23:05:15  caron
# upgrade to java2D, new ProjectionManager
#
# Revision 1.2  1998/12/14  17:10:44  russ
# Add comment for accumulating change histories.
#
*/







