//
// C++ Interface: earth_geometry
//
// Description: set of constants that defines the geometry of the earth
//
//
// Author: Nicolas PASCAL <nicolas.pascal@icare.univ-lille1.fr>, (C) 2012
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef EARTH_GEOMETRY_H
#define EARTH_GEOMETRY_H

#include "geometry_common.h"

/**
 * @namespace Earth geometry constants of the earth
 * Source : http://en.wikipedia.org/wiki/Earth_radius
 */
namespace Earth {
    /** gravity constant as defined by the WMO (m.s-2) */
    static const double g0 = 9.80665;
    /** if earth is considered as a perfect sphere, authalic radius ("equal aera") (m) */
    static const double authalic_radius = 6371007.2;

    /**
     * @enum Model define the earth model to use
     */
    typedef enum Model {
        SPHERICAL_MODEL = 0, // sphere
        WGS84_MODEL, // ellipsoid ; World Geodetic System 1984
    } Model;

    /**
     * @namespace WGS84 geometry constants of the earth, as defined by the World Geodetic System 1984 convention (ie the one used by the GPS system). In this case, earth is considered as an ellipsoid
     */
    namespace WGS84 {
        /** equatorial radius, semi-major axis (m) */
        const double a = 6378137.0;
        /** polare radius, semi-minor axis (m) */
        const double b = 6356752.314245;
        /** square of semi-major axis */
        const double a2 = a * a;
        /** 4th power of semi-major axis */
        const double a4 = a2 * a2;
        /** square of semi-minor axis */
        const double b2 = b * b;
        /** 4th power of semi-minor axis */
        const double b4 = b2 * b2;
        /** square of the eccentricity */
        const double e2 = 1. - (b2 / a2);
        /** eccentricity */
        const double e = sqrt (e2); // eccentricity
        const double ab =  a / b;
        const double e22 = (a2 - b2) / b2;
        const double Ime2 = 1. - e2;
        /**
        * @brief computes N the radius of curvature in the prime vertical at the given latitude, equal to a/sqrt(1 - e^2 sin^2 B)
        * @param theta the colatitude, in radians
        * @return radius of curvature in the prime vertical at the given latitude
        */
        double get_N(double theta);
    }
    /**
     * @brief return the distance to the center of the earth in meters
     * Based on the WGS84 conventions, if the earth is considered as :
     * - a sphere, it will have a radius of 6371007.2 m
     * - an ellipsoid, the equatorial radius will be 6378137.0 m and the polar 6356752.3142 m
     * @return the radius of the earth in meters, eventually at the given latitude
     */
    double get_radius ();
}

/**
 * @brief compute the inverse distance weighting interpolation ( IDW ) of a variable between 4 points on a sphere identified by their (lat, lon) coordinates and their value
 * @param lat latitude of the interpolation point P, in degrees and in the range [ -90, 90 ]
 * @param lon longitude of the interpolation point, in degrees and in the range ] -180, 180 ]
 * @param lat1 latitude of the 1st point
 * @param lon1 longitude of the 1st point
 * @param v1 value of the variable to interpolate on the 1st point
 * @param lat2 latitude of the 2nd point
 * @param lon2 longitude of the 2nd point
 * @param v2 value of the variable to interpolate on the 2nd point
 * @param lat3 latitude of the 3rd point
 * @param lon3 longitude of the 3rd point
 * @param v3 value of the variable to interpolate on the 3rd point
 * @param lat4 latitude of the 4th point
 * @param lon4 longitude of the 4th point
 * @param v4 value of the variable to interpolate on the 4th point
 * @return the interpolated value on the point (lat, lon)
 */
// inline const double IDW_interp ( double lat,  double lon,
//                                  double lat1, double lon1, double v1,
//                                  double lat2, double lon2, double v2,
//                                  double lat3, double lon3, double v3,
//                                  double lat4, double lon4, double v4
//                                ) {
//
//     // compute the haversine distance from P to P1, P2, P3 and P4
//     double d1 = Geocentric::get_haversine_dist ( lat, lon, lat1, lon1 );
//     double d2 = Geocentric::get_haversine_dist ( lat, lon, lat2, lon2 );
//     double d3 = Geocentric::get_haversine_dist ( lat, lon, lat3, lon3 );
//     double d4 = Geocentric::get_haversine_dist ( lat, lon, lat4, lon4 );
//
//     // compute the corresponding weights
//     double w1 = 1. / pow ( d1, 2. );
//     double w2 = 1. / pow ( d2, 2. );
//     double w3 = 1. / pow ( d3, 2. );
//     double w4 = 1. / pow ( d4, 2. );
//
//     // compute interpolation
//     return ( w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4 ) / ( w1 + w2 + w3 + w4 );
// };
//
// inline const double DW_interp (  double lat,  double lon,
//                                  double lat1, double lon1, double v1,
//                                  double lat2, double lon2, double v2,
//                                  double lat3, double lon3, double v3,
//                                  double lat4, double lon4, double v4
//                                ) {
//     // compute the haversine distance from P to P1, P2, P3 and P4
//     double d1 = Geocentric::get_haversine_dist ( lat, lon, lat1, lon1 );
//     double d2 = Geocentric::get_haversine_dist ( lat, lon, lat2, lon2 );
//     double d3 = Geocentric::get_haversine_dist ( lat, lon, lat3, lon3 );
//     double d4 = Geocentric::get_haversine_dist ( lat, lon, lat4, lon4 );
//
//     // compute the corresponding weights
//     double w1 = d1;
//     double w2 = d2;
//     double w3 = d3;
//     double w4 = d4;
//
//     // compute interpolation
//     return ( w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4 ) / ( w1 + w2 + w3 + w4 );
// };
#endif