(function() {
  var gju = this.gju = {};
  
  // Export the geojson object for **CommonJS**
  if (typeof module !== 'undefined' && module.exports) {
    module.exports = gju;
  }
  
  // adapted from http://www.kevlindev.com/gui/math/intersection/Intersection.js
  gju.lineStringsIntersect = function(l1, l2) {
  var intersects = [];
  for (var i = 0; i <= l1.coordinates.length - 2; ++i) {
    for (var j = 0; j <= l2.coordinates.length - 2; ++j) {      
    var a1 = {x: l1.coordinates[i][1], y: l1.coordinates[i][0]},
      a2 = {x: l1.coordinates[i+1][1], y: l1.coordinates[i+1][0]},
      b1 = {x: l2.coordinates[j][1], y: l2.coordinates[j][0]},
      b2 = {x: l2.coordinates[j+1][1], y: l2.coordinates[j+1][0]},
      ua_t = (b2.x - b1.x) * (a1.y - b1.y) - (b2.y - b1.y) * (a1.x - b1.x),
      ub_t = (a2.x - a1.x) * (a1.y - b1.y) - (a2.y - a1.y) * (a1.x - b1.x),
      u_b  = (b2.y - b1.y) * (a2.x - a1.x) - (b2.x - b1.x) * (a2.y - a1.y);
    if ( u_b != 0 ) {
      var ua = ua_t / u_b,
        ub = ub_t / u_b;
      if ( 0 <= ua && ua <= 1 && 0 <= ub && ub <= 1 ) {
      intersects.push({ 
        'type': 'Point',
        'coordinates': [a1.x + ua * (a2.x - a1.x), a1.y + ua * (a2.y - a1.y)]
      });
      }
    }
    }
  }
  if (intersects.length == 0) intersects = false;
    return intersects;
  }

  // adapted from http://jsfromhell.com/math/is-point-in-poly
  gju.pointInPolygon = function(point, polygon) {
  var x = point.coordinates[1],
    y = point.coordinates[0],
    poly = polygon.coordinates[0]; //TODO: support polygons with holes
  for (var c = false, i = -1, l = poly.length, j = l - 1; ++i < l; j = i) {
    var px = poly[i][1], py = poly[i][0],
        jx = poly[j][1], jy = poly[j][0];
    if (((py <= y && y < jy) || (jy <= y && y < py)) && (x < (jx - px) * (y - py) / (jy - py) + px)) {
      c = [point];
    }
  }
  return c;
  }
  
  gju.numberToRadius = function(number) {
    return number * Math.PI / 180;
  }

  gju.numberToDegree = function(number) {
    return number * 180 / Math.PI;
  }
  
  // written with help from @tautologe 
  gju.drawCircle = function(radiusInMeters, centerPoint) {
    var center = [centerPoint.coordinates[1], centerPoint.coordinates[0]],
      dist = (radiusInMeters / 1000) / 6371, // convert meters to radiant
      radCenter = [gju.numberToRadius(center[0]), gju.numberToRadius(center[1])],
      steps = 15, // 15 sided circle
      poly = [[center[0], center[1]]];
    for (var i = 0; i < steps + 1; i++) {
    	var brng = 2 * Math.PI * i / steps;  
    	var lat = Math.asin(Math.sin(radCenter[0]) * Math.cos(dist) + 
        Math.cos(radCenter[0]) * Math.sin(dist) * Math.cos(brng));
    	var lng = radCenter[1] + Math.atan2(Math.sin(brng) * Math.sin(dist) *
        Math.cos(radCenter[0]), 
        Math.cos(dist) - Math.sin(radCenter[0]) *
        Math.sin(lat));
      poly[i] = [];
      poly[i][1] = gju.numberToDegree(lat);
      poly[i][0] = gju.numberToDegree(lng);
    }
    return { "type": "Polygon", "coordinates": [poly] };
  }

  gju.rectangleCentroid = function(rectangle) {
    var bbox = rectangle.coordinates[0];
    var xmin = bbox[0][0], ymin = bbox[0][1], xmax = bbox[1][0], ymax = bbox[1][1];
    var xwidth = xmax - xmin;
    var ywidth = ymax - ymin;
    return { 'type': 'Point', 'coordinates': [xmin + xwidth/2, ymin + ywidth/2] };
  }
  
  // from http://www.movable-type.co.uk/scripts/latlong.html
  gju.pointDistance = function(pt1, pt2) {
    var lon1 = pt1.coordinates[0], lat1 = pt1.coordinates[1],
      lon2 = pt2.coordinates[0], lat2 = pt2.coordinates[1],
      dLat = gju.numberToRadius(lat2 - lat1),
      dLon = gju.numberToRadius(lon2 - lon1),
      a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.cos(gju.numberToRadius(lat1)) * Math.cos(gju.numberToRadius(lat2)) * 
        Math.sin(dLon/2) * Math.sin(dLon/2),
      c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    return (6371 * c) * 1000; // returns meters
  },

  // checks if geometry lies entirely within a circle
  // works with Point, LineString, Polygon
  gju.geometryWithinRadius = function(geometry, center, radius) {
    if (geometry.type == 'Point') {
      return gju.pointDistance(geometry, center) <= radius;
    } else if (geometry.type == 'LineString' || geometry.type == 'Polygon') {
      var point = {};
      var coordinates;
      if (geometry.type == 'Polygon') {
        // it's enough to check the exterior ring of the Polygon
        coordinates = geometry.coordinates[0];
      } else {
        coordinates = geometry.coordinates;
      }
      for (var i in coordinates) {
        point.coordinates = coordinates[i];
        if (gju.pointDistance(point, center) > radius) {
          return false;
        }
      }
    }
    return true;
  }
  
  // adapted from http://paulbourke.net/geometry/polyarea/javascript.txt
  gju.area = function(polygon) {
    var area = 0;
    // TODO: polygon holes at coordinates[1]
    var points = polygon.coordinates[0];
    var j = points.length - 1;
    var p1, p2;

    for (var i=0; i < points.length; j = i++) {
      var p1 = {x: points[i][1], y: points[i][0]};
      var p2 = {x: points[j][1], y: points[j][0]};
      area += p1.x * p2.y;
      area -= p1.y * p2.x;
    }

    area /= 2;
    return area;
  },

  // adapted from http://paulbourke.net/geometry/polyarea/javascript.txt
  gju.centroid = function(polygon) {
    var f, x = 0, y = 0;
    // TODO: polygon holes at coordinates[1]
    var points = polygon.coordinates[0];
    var j = points.length - 1;
    var p1, p2;

    for (var i=0; i < points.length; j = i++) {
      var p1 = {x: points[i][1], y: points[i][0]};
      var p2 = {x: points[j][1], y: points[j][0]};
      f = p1.x * p2.y - p2.x * p1.y;
      x += (p1.x + p2.x) * f;
      y += (p1.y + p2.y) * f;
    }

    f = gju.area(polygon) * 6;
    return { 'type': 'Point', 'coordinates': [y/f, x/f] };
  },

  gju.simplify = function (source, kink) {
    /* source[] array of geojson points */
    /* kink	in metres, kinks above this depth kept  */
    /* kink depth is the height of the triangle abc where a-b and b-c are two consecutive line segments */
    kink = kink || 20;
    source = source.map(function(o) { return {lng: o.coordinates[0], lat: o.coordinates[1]} });
  
    var	n_source, n_stack, n_dest, start, end, i, sig;  
    var dev_sqr, max_dev_sqr, band_sqr;
    var x12, y12, d12, x13, y13, d13, x23, y23, d23;
    var F = ((Math.PI / 180.0) * 0.5 );
    var index = new Array(); /* aray of indexes of source points to include in the reduced line */
    var sig_start = new Array(); /* indices of start & end of working section */
    var sig_end = new Array();	
  
    /* check for simple cases */
  
    if ( source.length < 3 ) return(source);  /* one or two points */
  
    /* more complex case. initialize stack */

  	n_source = source.length;
    band_sqr = kink * 360.0 / (2.0 * Math.PI * 6378137.0);	/* Now in degrees */
    band_sqr *= band_sqr;
    n_dest = 0;
    sig_start[0] = 0;
    sig_end[0] = n_source-1;
    n_stack = 1;

    /* while the stack is not empty  ... */
    while ( n_stack > 0 ){

      /* ... pop the top-most entries off the stacks */

      start = sig_start[n_stack-1];
      end = sig_end[n_stack-1];
      n_stack--;

      if ( (end - start) > 1 ){  /* any intermediate points ? */    

          /* ... yes, so find most deviant intermediate point to
             either side of line joining start & end points */                   

        x12 = (source[end].lng() - source[start].lng());
        y12 = (source[end].lat() - source[start].lat());
        if (Math.abs(x12) > 180.0) 
          x12 = 360.0 - Math.abs(x12);
        x12 *= Math.cos(F * (source[end].lat() + source[start].lat()));/* use avg lat to reduce lng */
        d12 = (x12*x12) + (y12*y12);

        for ( i = start + 1, sig = start, max_dev_sqr = -1.0; i < end; i++ ){                  

          x13 = (source[i].lng() - source[start].lng());
          y13 = (source[i].lat() - source[start].lat());
          if (Math.abs(x13) > 180.0) 
            x13 = 360.0 - Math.abs(x13);
          x13 *= Math.cos (F * (source[i].lat() + source[start].lat()));
          d13 = (x13*x13) + (y13*y13);

          x23 = (source[i].lng() - source[end].lng());
          y23 = (source[i].lat() - source[end].lat());
          if (Math.abs(x23) > 180.0) 
            x23 = 360.0 - Math.abs(x23);
          x23 *= Math.cos(F * (source[i].lat() + source[end].lat()));
          d23 = (x23*x23) + (y23*y23);

          if ( d13 >= ( d12 + d23 ) )
            dev_sqr = d23;
          else if ( d23 >= ( d12 + d13 ) )
            dev_sqr = d13;
          else
            dev_sqr = (x13 * y12 - y13 * x12) * (x13 * y12 - y13 * x12) / d12;// solve triangle

          if ( dev_sqr > max_dev_sqr  ){
            sig = i;
            max_dev_sqr = dev_sqr;
          }
        }

        if ( max_dev_sqr < band_sqr ){   /* is there a sig. intermediate point ? */
          /* ... no, so transfer current start point */
          index[n_dest] = start;
          n_dest++;
        }
        else{
          /* ... yes, so push two sub-sections on stack for further processing */
          n_stack++;
          sig_start[n_stack-1] = sig;
          sig_end[n_stack-1] = end;
          n_stack++;
          sig_start[n_stack-1] = start;
          sig_end[n_stack-1] = sig;
        }
      }
      else{
          /* ... no intermediate points, so transfer current start point */
          index[n_dest] = start;
          n_dest++;
      }
    }

    /* transfer last point */
    index[n_dest] = n_source-1;
    n_dest++;

    /* make return array */
    var r = new Array();
    for(var i=0; i < n_dest; i++)
      r.push(source[index[i]]);
    return r.map(function(o) { return {type: "Point", coordinates: [o.lng, o.lat]} });;
  }
})();