/*
 *  Hamlib Interface - locator and bearing conversion calls
 *  Copyright (c) 2001-2010 by Stephane Fillod
 *  Copyright (c) 2003 by Nate Bargmann
 *  Copyright (c) 2003 by Dave Hines
 *
 *
 *  Code to determine bearing and range was taken from the Great Circle,
 *  by S. R. Sampson, N5OWK.
 *  Ref: "Air Navigation", Air Force Manual 51-40, 1 February 1987
 *  Ref: "ARRL Satellite Experimenters Handbook", August 1990
 *
 *  Code to calculate distance and azimuth between two Maidenhead locators,
 *  taken from wwl, by IK0ZSN Mirko Caserta.
 *
 *  New bearing code added by N0NB was found at:
 *    http://williams.best.vwh.net/avform.htm#Crs
 *
 *
 *   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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 */


/**
 * \addtogroup utilities
 * @{
 */

/**
 * \file src/locator.c
 *
 * \brief QRA locator (Maidenhead grid square) and latitude/longitude bearing
 * conversion interface.
 *
 * \author Stephane Fillod
 * \author Nate Bargmann
 * \author Dave Hines
 * \author The Hamlib Group
 * \date 2000-2020
 */

/**
 * \page hamlib Hamlib general purpose API
 *
 * Hamlib function call interface for determining QRA locator (Maidenhead grid
 * square), bearing, and conversion between QRA locator and latitude/longitude
 * formats.
 *
 * \par Sources used in writing these routines
 *
 * \parblock
 * Code to determine bearing and range was taken from the Great Circle,
 * by Steven R. Sampson, N5OWK.<br />
 * Ref: "Air Navigation", Air Force Manual 51-40, 1 February 1987<br />
 * Ref: "ARRL Satellite Experimenters Handbook", August 1990
 *
 * Code to calculate distance and azimuth between two QRA locators, taken from
 * wwl, by IK0ZSN, Mirko Caserta.
 *
 * New bearing code added by N0NB was found at:
 * http://williams.best.vwh.net/avform.htm#Crs
 * \endparblock
 */


#include <hamlib/config.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>


#include <hamlib/rotator.h>


/** \brief Standard definition of a radian. */
#define RADIAN  (180.0 / M_PI)

/** \brief arc length for 1 degree in kilometers, i.e. 60 Nautical Miles */
#define ARC_IN_KM 111.2


/* The following is contributed by Dave Hines M1CXW
 *
 * begin dph
 */

/* At this time documenting a single static variable as in loc_char_range[]
 * below is not supported by Doxygen.  Hide this section until support exists
 * or a work-around becomes available.
 */
#ifndef DOC_HIDDEN

/**
 * \brief Constants used when converting between Maidenhead grid
 * locators and longitude/latitude values.
 *
 * \ref MAX_LOCATOR_PAIRS is the maximum number of locator character pairs to
 * convert. This number MUST NOT exceed the number of pairs of values in
 * loc_char_range[].  Setting \ref MAX_LOCATOR_PAIRS to 3 will convert the
 * currently defined 6 character locators.  A value of 4 will convert the
 * extended 8 character locators described in section 3L of "The IARU region 1
 * VHF managers handbook".  Values of 5 and 6 will extent the format even
 * more, to the longest definition I have seen for locators, see
 * http://www.btinternet.com/~g8yoa/geog/non-ra.html (currently a dead
 * link. -N0NB).  Be aware that there seems to be no universally accepted
 * standard for 10 & 12 character locators.
 *
 * The ranges of characters which will be accepted by locator2longlat(), and
 * generated by longlat2locator(), are specified by the \ref loc_char_range[]
 * array.  This array may be changed without requiring any other code changes.
 *
 * For the fifth pair to range from aa to xx use:
 * \code const static int loc_char_range[] = { 18, 10, 24, 10, 24, 10 };\endcode
 *
 * For the fifth pair to range from aa to yy use:
 * \code const static int loc_char_range[] = { 18, 10, 24, 10, 25, 10 };\endcode
 */
const static int loc_char_range[] = { 18, 10, 24, 10, 24, 10 };

#endif  /* !DOC_HIDDEN */

/** \def MAX_LOCATOR_PAIRS
 *
 * \brief Longest locator to process, e.g. AA00AA00AA00.
 *
 * Sets the limit locator2longlat() will convert and sets the maximum length
 * longlat2locator() will generate.  Each function properly handles any value
 * from `1` to `6` so \ref MAX_LOCATOR_PAIRS should be left at `6`.
 *
 * \def MIN_LOCATOR_PAIRS
 *
 * \brief Shortest locator to process, e.g. AA.
 *
 * Sets a floor on the shortest locator that should be handled.
 */
#define MAX_LOCATOR_PAIRS       6
#define MIN_LOCATOR_PAIRS       1

/* end dph */


/**
 * \brief Convert Degrees Minutes Seconds (DMS) notation to decimal degrees
 * (D.DDD) angle.
 *
 * \param degrees Degrees, whole degrees.
 * \param minutes Minutes, whole minutes.
 * \param seconds Seconds, decimal seconds.
 * \param sw South or West.
 *
 * Convert a Degrees Minutes Seconds (DMS) notation value to a decimal degrees
 * (D.DDD) angle value.
 *
 * \note For the parameters \a degrees >360, \a minutes > 60, and \a seconds >
 * 60.0 are allowed, but the resulting angle will not be normalized.
 *
 * When the variable \a sw is passed a value of 1, the returned decimal
 * degrees value will be negative (*South* or *West*).  When passed a value of 0
 * the returned decimal degrees value will be positive (*North* or *East*).
 *
 * \return The signed angle in decimal degrees (D.DDD).
 *
 * \sa dec2dms()
 */
double HAMLIB_API dms2dec(int degrees, int minutes, double seconds, int sw)
{
    double st;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    if (degrees < 0)
    {
        degrees = abs(degrees);
    }

    if (minutes < 0)
    {
        minutes = abs(minutes);
    }

    if (seconds < 0)
    {
        seconds = fabs(seconds);
    }

    st = (double)degrees + (double)minutes / 60. + seconds / 3600.;

    if (sw == 1)
    {
        return -st;
    }
    else
    {
        return st;
    }
}


/**
 * \brief Convert degrees decimal minutes (D M.MMM) notation to decimal
 * degrees (D.DDD) angle.
 *
 * \param degrees Degrees, whole degrees.
 * \param minutes Minutes, decimal minutes.
 * \param seconds Seconds, decimal seconds.
 * \param sw South or West.
 *
 * Convert a degrees decimal minutes (D M.MMM) notation common on many GPS
 * units to a decimal degrees (D.DDD) angle value.
 *
 * \note For the parameters \a degrees > 360, \a minutes > 60.0, \a seconds >
 * 60.0 are allowed, but the resulting angle will not be normalized.
 *
 * When the variable \a sw is passed a value of 1, the returned decimal
 * degrees value will be negative (*South* or *West*).  When passed a value of
 * 0 the returned decimal degrees value will be positive (*North* or *East*).
 *
 * \return The signed angle in decimal degrees (D.DDD).
 *
 * \sa dec2dmmm()
 */
double HAMLIB_API dmmm2dec(int degrees, double minutes, double seconds, int sw)
{
    double st;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    if (degrees < 0)
    {
        degrees = abs(degrees);
    }

    if (minutes < 0)
    {
        minutes = fabs(minutes);
    }

    st = (double)degrees + (minutes / 60) + (seconds / 3600);

    if (sw == 1)
    {
        return -st;
    }
    else
    {
        return st;
    }
}


/**
 * \brief Convert a decimal degrees (D.DDD) angle into Degrees Minutes
 * Seconds (DMS) notation.
 *
 * \param dec Decimal degrees (D.DDD).
 * \param degrees Pointer for the calculated whole Degrees.
 * \param minutes Pointer for the calculated whole Minutes.
 * \param seconds Pointer for the calculated decimal Seconds.
 * \param sw Pointer for the calculated SW (South/West) flag.
 *
 * Convert decimal degrees angle (D.DDD) into its Degree Minute Second (DMS)
 * notation.
 *
 * When \a dec < -180 or \a dec > 180, the angle will be normalized within
 * these limits and the sign set appropriately.
 *
 * Upon return, guarantees 0 >= \a degrees <= 180, 0 >= \a minutes < 60, and
 * 0.0 >= \a seconds < 60.0.
 *
 * When \a dec is < 0.0 \a sw will be set to 1.  When \a dec is >= 0.0 \a sw
 * will be set to 0.  This flag allows the application to determine whether
 * the DMS angle should be treated as negative (*South* or *West*).
 *
 * \return RIG_OK if the operation has been successful, otherwise a **negative
 * value** if an error occurred (in which case, cause is set appropriately).
 *
 * \retval RIG_OK The conversion was successful.
 * \retval RIG_EINVAL Either of the pointers are NULL.
 *
 * \sa dms2dec()
 */
int HAMLIB_API dec2dms(double dec,
                       int *degrees,
                       int *minutes,
                       double *seconds,
                       int *sw)
{
    int deg, min;
    double st;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    /* bail if NULL pointers passed */
    if (!degrees || !minutes || !seconds || !sw)
    {
        return -RIG_EINVAL;
    }

    /* reverse the sign if dec has a magnitude greater
     * than 180 and factor out multiples of 360.
     * e.g. when passed 270 st will be set to -90
     * and when passed -270 st will be set to 90.  If
     * passed 361 st will be set to 1, etc.  If passed
     * a value > -180 || < 180, value will be unchanged.
     */
    if (dec >= 0.0)
    {
        st = fmod(dec + 180, 360) - 180;
    }
    else
    {
        st = fmod(dec - 180, 360) + 180;
    }

    /* if after all of that st is negative, we want deg
     * to be negative as well except for 180 which we want
     * to be positive.
     */
    if (st < 0.0 && st != -180)
    {
        *sw = 1;
    }
    else
    {
        *sw = 0;
    }

    /* work on st as a positive value to remove a
     * bug introduced by the effect of floor() when
     * passed a negative value.  e.g. when passed
     * -96.8333 floor() returns -95!  Also avoids
     * a rounding error introduced on negative values.
     */
    st = fabs(st);

    deg = (int)floor(st);
    st  = 60. * (st - (double)deg);
    min = (int)floor(st);
    st  = 60. * (st - (double)min);

    *degrees = deg;
    *minutes = min;
    *seconds = st;

    return RIG_OK;
}


/**
 * \brief Convert a decimal degrees (D.DDD) angle into degrees decimal minutes
 * (D M.MMM) notation.
 *
 * \param dec Decimal degrees angle.
 * \param degrees Pointer for the calculated whole Degrees.
 * \param minutes Pointer for the calculated decimal Minutes.
 * \param sw Pointer for the calculated SW flag.
 *
 * Convert a decimal angle into its degree, decimal minute
 * notation common on many GPS units.
 *
 * When passed a value < -180 or > 180, the value will be normalized
 * within these limits and the sign set appropriately.
 *
 * Upon return dec2dmmm guarantees 0 >= \a degrees <= 180,
 * 0.0 >= \a minutes < 60.0.
 *
 * When \a dec is < 0.0 \a sw will be set to 1.  When \a dec is
 * >= 0.0 \a sw will be set to 0.  This flag allows the application
 * to determine whether the D M.MMM angle should be treated as negative
 * (south or west).
 *
 * \return RIG_OK if the operation has been successful, otherwise a **negative
 * value** if an error occurred (in which case, cause is set appropriately).
 *
 * \retval RIG_OK The conversion was successful.
 * \retval RIG_EINVAL Either of the pointers are NULL.
 *
 * \sa dmmm2dec()
 */
int HAMLIB_API dec2dmmm(double dec, int *degrees, double *minutes, int *sw)
{
    int r, min;
    double sec;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    /* bail if NULL pointers passed */
    if (!degrees || !minutes || !sw)
    {
        return -RIG_EINVAL;
    }

    r = dec2dms(dec, degrees, &min, &sec, sw);

    if (r != RIG_OK)
    {
        return r;
    }

    *minutes = (double)min + sec / 60;

    return RIG_OK;
}


/**
 * \brief Convert QRA locator (Maidenhead grid square) to Longitude/Latitude.
 *
 * \param longitude Pointer for the calculated Longitude.
 * \param latitude Pointer for the calculated Latitude.
 * \param locator The QRA locator--2 through 12 characters + nul string.
 *
 * Convert a QRA locator string to Longitude/Latitude in decimal degrees
 * (D.DDD).  The locator should be 2 through 12 chars long format.
 * \a locator2longlat is case insensitive, however it checks for locator
 * validity.
 *
 * Decimal long/lat is computed to center of grid square, i.e. given
 * `EM19` will return coordinates equivalent to the southwest corner
 * of `EM19mm`.
 *
 * \return RIG_OK if the operation has been successful, otherwise a **negative
 * value** if an error occurred (in which case, cause is set appropriately).
 *
 * \retval RIG_OK The conversion was successful.
 * \retval RIG_EINVAL The QRA locator exceeds RR99xx99xx99 or exceeds length
 * limit--currently 1 to 6 lon/lat pairs--or is otherwise malformed.
 *
 * \bug The fifth pair ranges from aa to xx, there is another convention
 *  that ranges from aa to yy.  At some point both conventions should be
 *  supported.
 *
 * \sa longlat2locator()
 */
/* begin dph */
int HAMLIB_API locator2longlat(double *longitude,
                               double *latitude,
                               const char *locator)
{
    int x_or_y, paircount;
    int locvalue, pair;
    double xy[2];

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    /* bail if NULL pointers passed */
    if (!longitude || !latitude)
    {
        return -RIG_EINVAL;
    }

    paircount = strlen(locator) / 2;

    /* verify paircount is within limits */
    if (paircount > MAX_LOCATOR_PAIRS)
    {
        paircount = MAX_LOCATOR_PAIRS;
    }
    else if (paircount < MIN_LOCATOR_PAIRS)
    {
        return -RIG_EINVAL;
    }

    /* For x(=longitude) and y(=latitude) */
    for (x_or_y = 0;  x_or_y < 2;  ++x_or_y)
    {
        double ordinate = -90.0;
        int divisions = 1;

        for (pair = 0;  pair < paircount;  ++pair)
        {
            locvalue = locator[pair * 2 + x_or_y];

            /* Value of digit or letter */
            locvalue -= (loc_char_range[pair] == 10) ? '0' :
                        (isupper(locvalue)) ? 'A' : 'a';

            /* Check range for non-letter/digit or out of range */
            if ((locvalue < 0) || (locvalue >= loc_char_range[pair]))
            {
                return -RIG_EINVAL;
            }

            divisions *= loc_char_range[pair];
            ordinate += locvalue * 180.0 / divisions;
        }

        /* Center ordinate in the Maidenhead "square" or "subsquare" */
        ordinate += 90.0 / divisions;

        xy[x_or_y] = ordinate;
    }

    *longitude = xy[0] * 2.0;
    *latitude = xy[1];

    return RIG_OK;
}
/* end dph */


/**
 * \brief Convert longitude/latitude to QRA locator (Maidenhead grid square).
 *
 * \param longitude Longitude, decimal degrees.
 * \param latitude Latitude, decimal degrees.
 * \param locator Pointer for the QRA Locator.
 * \param pair_count Requested precision expressed as lon/lat pairs in the
 * returned QRA locator string.
 *
 * Convert longitude/latitude given in decimal degrees (D.DDD) to a QRA
 * locator (Maidenhead grid square).  \a locator must point to an array length
 * that is at least \a pair_count * 2 char + '\\0'.
 *
 * \return RIG_OK if the operation has been successful, otherwise a **negative
 * value** if an error occurred (in which case, cause is set appropriately).
 *
 * \retval RIG_OK The conversion was successful.
 * \retval RIG_EINVAL if \a locator is NULL or \a pair_count exceeds length
 * limit.  Currently 1 to 6 lon/lat pairs.
 *
 * \bug \a locator is not tested for overflow.
 * \bug The fifth pair ranges from aa to yy, there is another convention
 * that ranges from aa to xx.  At some point both conventions should be
 * supported.
 *
 * \sa locator2longlat()
 */
/* begin dph */
int HAMLIB_API longlat2locator(double longitude,
                               double latitude,
                               char *locator,
                               int pair_count)
{
    int x_or_y, pair, locvalue;
    double square_size;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    if (!locator)
    {
        return -RIG_EINVAL;
    }

    if (pair_count < MIN_LOCATOR_PAIRS || pair_count > MAX_LOCATOR_PAIRS)
    {
        return -RIG_EINVAL;
    }

    for (x_or_y = 0;  x_or_y < 2;  ++x_or_y)
    {
        double ordinate = (x_or_y == 0) ? longitude / 2.0 : latitude;
        int divisions = 1;

        /* The 1e-6 here guards against floating point rounding errors */
        ordinate = fmod(ordinate + 270.000001, 180.0);

        for (pair = 0;  pair < pair_count;  ++pair)
        {
            divisions *= loc_char_range[pair];
            square_size = 180.0 / divisions;

            locvalue = (int)(ordinate / square_size);
            ordinate -= square_size * locvalue;
            locvalue += (loc_char_range[pair] == 10) ? '0' : 'A';
            locator[pair * 2 + x_or_y] = locvalue;
        }
    }

    locator[pair_count * 2] = '\0';

    return RIG_OK;
}
/* end dph */


/**
 * \brief Calculate the distance and bearing between two points.
 *
 * \param lon1 The local Longitude, decimal degrees.
 * \param lat1 The local Latitude, decimal degrees,
 * \param lon2 The remote Longitude, decimal degrees.
 * \param lat2 The remote Latitude, decimal degrees.
 * \param distance Pointer for the distance, km.
 * \param azimuth Pointer for the bearing, decimal degrees.
 *
 * Calculate the distance and bearing (QRB) between \a lon1, \a lat1 and
 * \a lon2, \a lat2.
 *
 * This version will calculate the QRB to a precision sufficient for 12
 * character locators.  Antipodal points, which are easily calculated, are
 * considered equidistant and the bearing is simply resolved to be true north,
 * e.g. \a azimuth = 0.0.
 *
 * \return RIG_OK if the operation has been successful, otherwise a **negative
 * value** if an error occurred (in which case, cause is set appropriately).
 *
 * \retval RIG_OK The calculations were successful.
 * \retval RIG_EINVAL If a NULL pointer passed or \a lat and \a lon values
 * exceed -90 to 90 or -180 to 180.
 *
 * \sa distance_long_path(), azimuth_long_path()
 */
int HAMLIB_API qrb(double lon1,
                   double lat1,
                   double lon2,
                   double lat2,
                   double *distance,
                   double *azimuth)
{
    double delta_long, tmp, arc, az;

    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    /* bail if NULL pointers passed */
    if (!distance || !azimuth)
    {
        return -RIG_EINVAL;
    }

    if ((lat1 > 90.0 || lat1 < -90.0) || (lat2 > 90.0 || lat2 < -90.0))
    {
        return -RIG_EINVAL;
    }

    if ((lon1 > 180.0 || lon1 < -180.0) || (lon2 > 180.0 || lon2 < -180.0))
    {
        return -RIG_EINVAL;
    }

    /* Prevent ACOS() Domain Error */
    if (lat1 == 90.0)
    {
        lat1 = 89.999999999;
    }
    else if (lat1 == -90.0)
    {
        lat1 = -89.999999999;
    }

    if (lat2 == 90.0)
    {
        lat2 = 89.999999999;
    }
    else if (lat2 == -90.0)
    {
        lat2 = -89.999999999;
    }

    /* Convert variables to Radians */
    lat1 /= RADIAN;
    lon1 /= RADIAN;
    lat2 /= RADIAN;
    lon2 /= RADIAN;

    delta_long = lon2 - lon1;

    tmp = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(delta_long);

    if (tmp > .999999999999999)
    {
        /* Station points coincide, use an Omni! */
        *distance = 0.0;
        *azimuth = 0.0;
        return RIG_OK;
    }

    if (tmp < -.999999)
    {
        /*
         * points are antipodal, it's straight down.
         * Station is equal distance in all Azimuths.
         * So take 180 Degrees of arc times 60 nm,
         * and you get 10800 nm, or whatever units...
         */
        *distance = 180.0 * ARC_IN_KM;
        *azimuth = 0.0;
        return RIG_OK;
    }

    arc = acos(tmp);

    /*
     * One degree of arc is 60 Nautical miles
     * at the surface of the earth, 111.2 km, or 69.1 sm
     * This method is easier than the one in the handbook
     */
    *distance = ARC_IN_KM * RADIAN * arc;

    /* Short Path */
    /* Change to azimuth computation by Dave Freese, W1HKJ */
    az = RADIAN * atan2(sin(lon2 - lon1) * cos(lat2),
                        (cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1)));

    az = fmod(360.0 + az, 360.0);

    if (az < 0.0)
    {
        az += 360.0;
    }
    else if (az >= 360.0)
    {
        az -= 360.0;
    }

    *azimuth = floor(az + 0.5);

    return RIG_OK;
}


/**
 * \brief Calculate the long path distance between two points.
 *
 * \param distance The shortpath distance in kilometers.
 *
 * Calculate the long path (opposite bearing of the short path) of a given
 * distance.
 *
 * \return The distance in kilometers for the opposite path.
 *
 * \sa qrb()
 */
double HAMLIB_API distance_long_path(double distance)
{
    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    return (ARC_IN_KM * 360.0) - distance;
}


/**
 * \brief Calculate the long path bearing between two points.
 *
 * \param azimuth The shortpath bearing--0.0 to 360.0 degrees.
 *
 * Calculate the long path (opposite of the short path) of a given bearing.
 *
 * \return the azimuth in decimal degrees for the opposite path or RIG_EINVAL
 * (negated value) upon input error (outside the range of 0.0 to 360.0).
 *
 * \sa qrb()
 */
double HAMLIB_API azimuth_long_path(double azimuth)
{
    rot_debug(RIG_DEBUG_VERBOSE, "%s called\n", __func__);

    if (azimuth == 0.0 || azimuth == 360.0)
    {
        return 180.0;
    }
    else if (azimuth > 0.0 && azimuth < 180.0)
    {
        return 180.0 + azimuth;
    }
    else if (azimuth == 180.0)
    {
        return 0.0;
    }
    else if (azimuth > 180.0 && azimuth < 360.0)
    {
        return (180.0 - azimuth) * -1.0;
    }
    else
    {
        return -RIG_EINVAL;
    }
}

/*! @} */