Location_lib with functions for area search

process/Location_lib/Location_lib.aod

+<?xml version="1.0" encoding="UTF-8"?>
+<process xmlns="http://www.adito.de/2018/ao/Model" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" VERSION="1.2.1" xsi:schemaLocation="http://www.adito.de/2018/ao/Model adito://models/xsd/process/1.2.1">
+  <name>Location_lib</name>
+  <majorModelMode>DISTRIBUTED</majorModelMode>
+  <process>%aditoprj%/process/Location_lib/process.js</process>
+  <variants>
+    <element>LIBRARY</element>
+  </variants>
+</process>

process/Location_lib/process.js

+import("Util_lib");
+import("system.eMath");
+
+function GeoLocationUtils () {}
+
+/**
+ * average earth radius in kilometers
+ */
+GeoLocationUtils.EARTH_RADIUS_AVG = 6371.001;
+
+/**
+ * earth radius at the equator in kilometers
+ */
+GeoLocationUtils.EARTH_RADIUS_EQUATOR = 6378.137;
+
+/**
+ * earth radius at the poles in kilometers
+ */
+GeoLocationUtils.EARTH_RADIUS_POLE = 6356.752;
+
+/**
+ * Calculates the distance between two given locations using the haversine formula.
+ * This method is reasonably accurate for most applications (in most cases the error is under 0.3%).
+ * 
+ * @param {String|Number} pLatA latitude of the first location
+ * @param {String|Number} pLonA longitude of the first location
+ * @param {String|Number} pLatB latitude of the second location
+ * @param {String|Number} pLonB longitude of the second location
+ * @return {Number} the distance in kilometers
+ */
+GeoLocationUtils.distanceHaversine = function (pLatA, pLonA, pLatB, pLonB)
+{
+    pLatA = GeoLocationUtils.degreeToRadians(pLatA);
+    pLonA = GeoLocationUtils.degreeToRadians(pLonA);
+    pLatB = GeoLocationUtils.degreeToRadians(pLatB);
+    pLonB = GeoLocationUtils.degreeToRadians(pLonB);
+
+    var R = GeoLocationUtils.EARTH_RADIUS_AVG; //average Earth radius in km https://en.wikipedia.org/wiki/Earth_radius
+
+    var deltaLat = eMath.subDec(pLatA, eMath.absDec(pLatB)); //lat can be negative
+    var deltaLon = eMath.subDec(pLonA, pLonB);
+
+    var sinLat = GeoLocationUtils.haversine(deltaLat);
+    var sinLon = GeoLocationUtils.haversine(deltaLon);
+    
+    return R * 2 * Math.asin(Math.sqrt(Math.pow(sinLat, 2) + Math.cos(pLatA) * Math.cos(eMath.absDec(pLatB)) * Math.pow(sinLon, 2)));
+}
+
+/**
+ * JavaScript function to calculate the geodetic distance between two points specified by latitude/longitude using the Vincenty inverse formula for ellipsoids.
+ * This is more accurate than the haversine formula but also more complex and thus theoretically not as fast.
+ *
+ * Taken from http://movable-type.co.uk/scripts/latlong-vincenty.html
+ * 
+ * @param   {Number} latA: first point in decimal degrees
+ * @param   {Number} lonA: first point in decimal degrees
+ * @param   {Number} latB: second point in decimal degrees
+ * @param   {Number} lonB: second point in decimal degrees
+ * @returns {Number} distance in kilometres between points
+ */
+GeoLocationUtils.distanceVincenty = function (latA, lonA, latB, lonB) 
+{
+    let phi1 = GeoLocationUtils.degreeToRadians(latA), lambda1 = GeoLocationUtils.degreeToRadians(lonA);
+    let phi2 = GeoLocationUtils.degreeToRadians(latB), lambda2 = GeoLocationUtils.degreeToRadians(lonB);
+
+    let a_ = 6378137;
+    let b = 6356752.314245;
+    let f = 1/298.257223563;
+    let epsilon = Math.pow(2, -52);
+
+    let L = lambda2 - lambda1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanphi.
+    let tanU1 = (1-f) * Math.tan(phi1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
+    let tanU2 = (1-f) * Math.tan(phi2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2;
+
+    let antipodal = Math.abs(L) > Math.PI/2 || Math.abs(phi2-phi1) > Math.PI/2;
+
+    let maxIterations = 100;
+    
+    let lambda = L, sinlambda = null, coslambda = null; // lambda = difference in longitude on an auxiliary sphere
+    let sigma = antipodal ? Math.PI : 0, sinsigma = 0, cossigma = antipodal ? -1 : 1, sinSqsigma = null; // sigma = angular distance P₁ P₂ on the sphere
+    let cos2sigmaM = 1;                      // sigmaM = angular distance on the sphere from the equator to the midpoint of the line
+    let sinalpha = null, cosSqalpha = 1;         // alpha = azimuth of the geodesic at the equator
+    let C = null;
+
+    let lambdaP = null, iterations = 0;
+    do {
+        sinlambda = Math.sin(lambda);
+        coslambda = Math.cos(lambda);
+        sinSqsigma = (cosU2*sinlambda) * (cosU2*sinlambda) + (cosU1*sinU2-sinU1*cosU2*coslambda) * (cosU1*sinU2-sinU1*cosU2*coslambda);
+        
+        if (Math.abs(sinSqsigma) < epsilon) 
+            break;  // co-incident/antipodal points (falls back on lambda/sigma = L)
+        
+        sinsigma = Math.sqrt(sinSqsigma);
+        cossigma = sinU1*sinU2 + cosU1*cosU2*coslambda;
+        sigma = Math.atan2(sinsigma, cossigma);
+        sinalpha = cosU1 * cosU2 * sinlambda / sinsigma;
+        cosSqalpha = 1 - sinalpha*sinalpha;
+        cos2sigmaM = (cosSqalpha != 0) ? (cossigma - 2*sinU1*sinU2/cosSqalpha) : 0; // on equatorial line cos²alpha = 0 (§6)
+        C = f/16*cosSqalpha*(4+f*(4-3*cosSqalpha));
+        lambdaP = lambda;
+        lambda = L + (1-C) * f * sinalpha * (sigma + C*sinsigma*(cos2sigmaM+C*cossigma*(-1+2*cos2sigmaM*cos2sigmaM)));
+        
+    } while (Math.abs(lambda-lambdaP) > 1e-12 && ++iterations<maxIterations);
+
+    if (iterations >= maxIterations)
+        return NaN;
+
+    let uSq = cosSqalpha * (a_*a_ - b*b) / (b*b);
+    let A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
+    let B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
+    let deltasigma = B*sinsigma*(cos2sigmaM+B/4*(cossigma*(-1+2*cos2sigmaM*cos2sigmaM)-
+        B/6*cos2sigmaM*(-3+4*sinsigma*sinsigma)*(-3+4*cos2sigmaM*cos2sigmaM)));
+
+    let s = b*A*(sigma-deltasigma); // s = length of the geodesic
+
+    return s / 1000; //convert metres into kilometres
+}
+
+/** 
+ * haversine function
+ * 
+ * @param {String|Number} pDelta req delta-value in radians
+ * @return {Number} hav(pDelta)
+ */ 
+GeoLocationUtils.haversine = function (pDelta)
+{
+    return Math.sin(eMath.divDec(pDelta, 2));
+}
+
+/**
+ * converts degrees to radians
+ * 
+ * @param {String|Number} pAngle angle in degrees
+ * @return {Number} angle in radians
+ */
+GeoLocationUtils.degreeToRadians = function (pAngle)
+{
+    return pAngle * Math.PI / 180;
+}
+
+function AreaSearchUtils () {}
+
+/**
+ * WIP
+ * @ignore
+ */
+AreaSearchUtils.findOrganisationsInArea = function (pLatitude, pLongitude, pMaxDistance)
+{
+    var locationData = [];
+    var locations = [];
+    
+    for (let i = 0, l = locationData.length; i < l; i++)
+    {
+        let distance = GeoLocationUtils.distanceVincenty(pLatitude, pLongitude, locationData[i][0], locationData[i][1]);
+        if (distace <= pMaxDistance)
+            locations.push(locationData[i]);
+    }
+    
+    return locations.sort(distanceSortFn);
+    
+    function distanceSortFn (a, b) 
+    { 
+        return a[3] - b[3]; 
+    }
+}
\ No newline at end of file