import { coordAll, featureEach } from '@turf/meta';
import { getCoords } from '@turf/invariant';
import { featureCollection, isObject, isNumber } from '@turf/helpers';
import centerMean from '@turf/center-mean';
import pointsWithinPolygon from '@turf/points-within-polygon';
import ellipse from '@turf/ellipse';

/**
 * Takes a {@link FeatureCollection} and returns a standard deviational ellipse,
 * also known as a “directional distribution.” The standard deviational ellipse
 * aims to show the direction and the distribution of a dataset by drawing
 * an ellipse that contains about one standard deviation’s worth (~ 70%) of the
 * data. 
 * 
 * This module was created in consultation with the following articles:
 *
 * • Robert S. Yuill, “The Standard Deviational Ellipse; An Updated Tool for 
 * Spatial Description,” _Geografiska Annaler_ 53, no. 1 (1971): 28–39, 
 * doi:{@link https://doi.org/10.2307/490885|10.2307/490885}.
 *
 * • Paul Hanly Furfey, “A Note on Lefever’s “Standard Deviational Ellipse,” 
 * _American Journal of Sociology_ 33, no. 1 (1927): 94—98, 
 * doi:{@link https://doi.org/10.1086/214336|10.1086/214336}.
 *
 * It mirrors the functionality of {@link http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-statistics-toolbox/directional-distribution.htm|Directional Distribution} in ArcGIS and the {@link http://arken.nmbu.no/~havatv/gis/qgisplugins/SDEllipse/|QGIS Standard Deviational Ellipse Plugin}.
 *
 * @name standardDeviationalEllipse
 * @param {Geometry|FeatureCollection<Point>} points GeoJSON points
 * @param {Object} [options={}] Optional parameters
 * @param {number} [options.weight] the property name used to weight the center
 * @param {number} [options.steps=64] number of steps for the polygon
 * @param {Object} [options.properties={}] properties to pass to the resulting ellipse
 * @returns {Feature<Polygon>} an elliptical Polygon that includes approximately 1 SD of the dataset within it.
 * @example
 *
 * var bbox = [-74, 40.72, -73.98, 40.74];
 * var points = turf.randomPoint(400, {bbox: bbox});
 * var sdEllipse = turf.standardDeviationalEllipse(points, options);
 * 
 * //addToMap
 * var addToMap = [points, sdEllipse];
 *
 */
function standardDeviationalEllipse(points, options) {
    // Optional params
    options = options || {};
    if (!isObject(options)) throw new Error('options is invalid');
    var steps = options.steps || 64;
    var weightTerm = options.weight;
    var properties = options.properties || {};

    // Validation:
    if (!isNumber(steps)) throw new Error('steps must be a number');
    if (!isObject(properties)) throw new Error('properties must be a number');

    // Calculate mean center:
    var meanCenter = centerMean(points, {weight: weightTerm});

    // Calculate angle of rotation:
    // [X, Y] = mean center of all [x, y].
    // theta = arctan( (A + B) / C )
    // A = sum((x - X)^2) - sum((y - Y)^2)
    // B = sqrt(A^2 + 4(sum((x - X)(y - Y))^2))
    // C = 2(sum((x - X)(y - Y)))

    var xDeviationSquaredSum = 0;
    var yDeviationSquaredSum = 0;
    var xyDeviationSum = 0;

    featureEach(points, function(point){
      var weight = point.properties[weightTerm] || 1;
      var deviation = getDeviations(getCoords(point), getCoords(meanCenter)); 
      xDeviationSquaredSum += Math.pow(deviation.x, 2) * weight;
      yDeviationSquaredSum += Math.pow(deviation.y, 2) * weight;
      xyDeviationSum += deviation.x * deviation.y * weight;
    });

    var bigA = xDeviationSquaredSum - yDeviationSquaredSum;
    var bigB = Math.sqrt(Math.pow(bigA, 2) + 4 * Math.pow(xyDeviationSum, 2));
    var bigC = 2 * xyDeviationSum;
    var theta = Math.atan((bigA + bigB) / bigC);
    var thetaDeg = theta * 180 / Math.PI;

    // Calculate axes:
    // sigmaX = sqrt((1 / n - 2) * sum((((x - X) * cos(theta)) - ((y - Y) * sin(theta)))^2))
    // sigmaY = sqrt((1 / n - 2) * sum((((x - X) * sin(theta)) - ((y - Y) * cos(theta)))^2))
    var sigmaXsum = 0;
    var sigmaYsum = 0;
    var weightsum = 0;
    featureEach(points, function(point){
      var weight = point.properties[weightTerm] || 1;
      var deviation = getDeviations(getCoords(point), getCoords(meanCenter)); 
      sigmaXsum += Math.pow((deviation.x * Math.cos(theta)) - (deviation.y * Math.sin(theta)), 2) * weight;
      sigmaYsum += Math.pow((deviation.x * Math.sin(theta)) - (deviation.y * Math.cos(theta)), 2) * weight;
      weightsum += weight;
    });

    var sigmaX = Math.sqrt(2 * sigmaXsum / weightsum);
    var sigmaY = Math.sqrt(2 * sigmaYsum / weightsum);

    var theEllipse = ellipse(meanCenter, sigmaX, sigmaY, {units: 'degrees', angle: thetaDeg, steps: steps, properties: properties});
    // var eccentricity = (Math.sqrt(Math.pow(sigmaX, 2) - Math.pow(sigmaY, 2))) / sigmaX;
    var pointsWithinEllipse = pointsWithinPolygon(points, featureCollection([theEllipse]));
    var standardDeviationalEllipseProperties = {
        meanCenterCoordinates: getCoords(meanCenter),
        semiMajorAxis: sigmaX,
        semiMinorAxis: sigmaY,
        angle: thetaDeg,
        // eccentricity: eccentricity,
        pctWithinEllipse: 100 * coordAll(pointsWithinEllipse).length / points.features.length
    };
    theEllipse.properties.standardDeviationalEllipse = standardDeviationalEllipseProperties;

    return theEllipse;
};

/**
 * Get x_i - X and y_i - Y
 *
 * @private
 * @param {array} coordinates Array of [x_i, y_i]
 * @param {array} center Array of [X, Y]
 * @returns {object} { x: n, y: m }
 */
function getDeviations(coordinates, center) {
    return {
      x: coordinates[0] - center[0], 
      y: coordinates[1] - center[1]
    };
}


export default standardDeviationalEllipse;

@@ -4330,33 +4330,33 @@
 ath.sin(theta)) 
--
++
  (deviation.y * 

Correct axis arithmetic error.