#ifndef GERMANAIRLINESVA_UTILITIES_GEODATA_H
#define GERMANAIRLINESVA_UTILITIES_GEODATA_H

#define M_PI 3.14159265358979323846
#define BUFSIZE 1024
#define MD5LEN 16
#define EARTH_M 6371000

#include <cmath>

namespace germanairlinesva
{
namespace utilities
{
  namespace geodata
  {
#pragma pack(push)
#pragma pack(1)
    struct point {
        double latitude;
        double longitude;
    };
#pragma pack(pop)

#pragma pack(push)
#pragma pack(1)
    struct box {
        struct point topLeft;
        struct point topRight;
        struct point bottomRight;
        struct point bottomLeft;
    };
#pragma pack(pop)

    static inline double toMetre(double value) { return value * 0.3048; }

    static inline double toFeet(double value) { return value * 3.280839895; }

    static inline double toDegrees(double value) { return value * 180 / M_PI; }

    static inline double toRadians(double value) { return value * M_PI / 180; }

    static inline double normalizeD(double value)
    {
      return fmod(value + 360, 360);
    }

    static inline double normalizeR(double value)
    {
      return value < 0 ? 2 * M_PI + value : value;
    }

    // Input and Output in degrees
    static inline double bearingDD(double fromLatitude,
                                   double fromLongitude,
                                   double toLatitude,
                                   double toLongitude)
    {
      double y = sin(toRadians(toLongitude) - toRadians(fromLongitude)) *
                 cos(toRadians(toLatitude));
      double x = cos(toRadians(fromLatitude)) * sin(toRadians(toLatitude)) -
                 sin(toRadians(fromLatitude)) * cos(toRadians(toLatitude)) *
                     cos(toRadians(toLongitude) - toRadians(fromLongitude));

      return normalizeD(toDegrees(atan2(y, x)));
    }

    // Input in degrees, Output in radians
    static inline double bearingDR(double fromLatitude,
                                   double fromLongitude,
                                   double toLatitude,
                                   double toLongitude)
    {
      double y = sin(toRadians(toLongitude) - toRadians(fromLongitude)) *
                 cos(toRadians(toLatitude));
      double x = cos(toRadians(fromLatitude)) * sin(toRadians(toLatitude)) -
                 sin(toRadians(fromLatitude)) * cos(toRadians(toLatitude)) *
                     cos(toRadians(toLongitude) - toRadians(fromLongitude));

      return normalizeR(atan2(y, x));
    }

    // Input in degrees, Output in metres
    static inline double distanceEarthD(double fromLatitude,
                                        double fromLongitude,
                                        double toLatitude,
                                        double toLongitude)
    {
      double lat1r, lon1r, lat2r, lon2r, u, v;
      lat1r = toRadians(fromLatitude);
      lon1r = toRadians(fromLongitude);
      lat2r = toRadians(toLatitude);
      lon2r = toRadians(toLongitude);
      u = sin((lat2r - lat1r) / 2);
      v = sin((lon2r - lon1r) / 2);
      return 2.0 * EARTH_M *
             asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
    }

    // Input in degrees, Output in metres
    static inline double distanceEarthP(point from, point to)
    {
      double lat1r, lon1r, lat2r, lon2r, u, v;
      lat1r = toRadians(from.latitude);
      lon1r = toRadians(from.longitude);
      lat2r = toRadians(to.latitude);
      lon2r = toRadians(to.longitude);
      u = sin((lat2r - lat1r) / 2);
      v = sin((lon2r - lon1r) / 2);
      return 2.0 * EARTH_M *
             asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
    }

    // Input and Output in degrees
    static inline struct point calculatePointDD(struct point coordinates,
                                                double bearing,
                                                double distance)
    {
      double r_latitude = toRadians(coordinates.latitude);
      double r_longitude = toRadians(coordinates.longitude);
      double r_bearing = toRadians(bearing);

      double newLatitude =
          std::asin(std::sin(r_latitude) * std::cos(distance / EARTH_M) +
                    std::cos(r_latitude) * std::sin(distance / EARTH_M) *
                        std::cos(r_bearing));
      double newLongitude =
          r_longitude +
          std::atan2(std::sin(r_bearing) * std::sin(distance / EARTH_M) *
                         std::cos(r_latitude),
                     std::cos(distance / EARTH_M) -
                         std::sin(r_latitude) * std::sin(newLatitude));

      return {toDegrees(newLatitude), toDegrees(newLongitude)};
    }

    // Input in degrees, Output in radians
    static inline struct point calculatePointDR(struct point coordinates,
                                                double bearing,
                                                double distance)
    {
      double r_latitude = toRadians(coordinates.latitude);
      double r_longitude = toRadians(coordinates.longitude);
      double r_bearing = toRadians(bearing);

      double newLatitude =
          std::asin(std::sin(r_latitude) * std::cos(distance / EARTH_M) +
                    std::cos(r_latitude) * std::sin(distance / EARTH_M) *
                        std::cos(r_bearing));
      double newLongitude =
          r_longitude +
          std::atan2(std::sin(r_bearing) * std::sin(distance / EARTH_M) *
                         std::cos(r_latitude),
                     std::cos(distance / EARTH_M) -
                         std::sin(r_latitude) * std::sin(newLatitude));

      return {newLatitude, newLongitude};
    }

    // LAT/LON in radians, BEARING/Output in degrees
    static inline struct point calculatePointRD(struct point coordinates,
                                                double bearing,
                                                double distance)
    {
      double r_bearing = toRadians(bearing);

      double newLatitude = std::asin(
          std::sin(coordinates.latitude) * std::cos(distance / EARTH_M) +
          std::cos(coordinates.latitude) * std::sin(distance / EARTH_M) *
              std::cos(r_bearing));
      double newLongitude =
          coordinates.longitude +
          std::atan2(std::sin(r_bearing) * std::sin(distance / EARTH_M) *
                         std::cos(coordinates.latitude),
                     std::cos(distance / EARTH_M) -
                         std::sin(coordinates.latitude) *
                             std::sin(newLatitude));

      return {toDegrees(newLatitude), toDegrees(newLongitude)};
    }

    // LAT/LON/Output in radians, BEARING in degrees
    static inline struct point calculatePointRR(struct point coordinates,
                                                double bearing,
                                                double distance)
    {
      double r_bearing = toRadians(bearing);

      double newLatitude = std::asin(
          std::sin(coordinates.latitude) * std::cos(distance / EARTH_M) +
          std::cos(coordinates.latitude) * std::sin(distance / EARTH_M) *
              std::cos(r_bearing));
      double newLongitude =
          coordinates.longitude +
          std::atan2(std::sin(r_bearing) * std::sin(distance / EARTH_M) *
                         std::cos(coordinates.latitude),
                     std::cos(distance / EARTH_M) -
                         std::sin(coordinates.latitude) *
                             std::sin(newLatitude));

      return {newLatitude, newLongitude};
    }

    // Input in degrees, calculate from center
    static inline box calculateBoxDD(struct point center,
                                     double bearing,
                                     double length,
                                     double width)
    {
      struct point primaryCenter =
          calculatePointDR(center, bearing, -length / 2);
      struct point secondaryCenter =
          calculatePointDR(center, bearing, length / 2);
      double offsetHeadingNorth =
          bearing - 90 > 0 ? bearing - 90 : bearing - 90 + 360;
      struct point primaryTop =
          calculatePointRD(primaryCenter, offsetHeadingNorth, width / 2);
      struct point primaryBottom =
          calculatePointRD(primaryCenter, offsetHeadingNorth, -width / 2);
      struct point secondaryTop =
          calculatePointRD(secondaryCenter, offsetHeadingNorth, width / 2);
      struct point secondaryBottom =
          calculatePointRD(secondaryCenter, offsetHeadingNorth, -width / 2);

      return {primaryTop, secondaryTop, secondaryBottom, primaryBottom};
    }

    // Input in degrees, calculate from start end
    static inline box calculateBoxDD(struct point start,
                                     struct point end,
                                     double bearing,
                                     double width)
    {
      double offsetHeadingNorth =
          bearing - 90 > 0 ? bearing - 90 : bearing - 90 + 360;
      struct point primaryTop =
          calculatePointDD(start, offsetHeadingNorth, width / 2);
      struct point primaryBottom =
          calculatePointDD(start, offsetHeadingNorth, -width / 2);
      struct point secondaryTop =
          calculatePointDD(end, offsetHeadingNorth, width / 2);
      struct point secondaryBottom =
          calculatePointDD(end, offsetHeadingNorth, -width / 2);

      return {primaryTop, secondaryTop, secondaryBottom, primaryBottom};
    }
  } // namespace geodata
} // namespace utilities
} // namespace germanairlinesva

#endif