/** 
Various utility methods 
*/

public class Utility {
	
    final static float EARTH_RADIUS_METERS = 6371000;

    /**
     * Calculate geodetic distance between two points specified by latitude/longitude using 
     * Vincenty inverse formula for ellipsoids; adapted from JavaScript at 
     * http://www.movable-type.co.uk/scripts/latlong-vincenty.html
     *
     * @param   {Number} lat1, lon1: first point in decimal degrees
     * @param   {Number} lat2, lon2: second point in decimal degrees
     * @returns (Number} distance in metres between points
     */
    public static double geodeticDistance(double lat1Deg, double lon1Deg, double lat2Deg, double lon2Deg) {
	double a = 6378137, b = 6356752.314245,  f = 1/298.257223563;  // WGS-84 ellipsoid params
	double L = Math.toRadians(lon2Deg-lon1Deg);
	double U1 = Math.atan((1-f) * Math.tan(Math.toRadians(lat1Deg)));
	double U2 = Math.atan((1-f) * Math.tan(Math.toRadians(lat2Deg)));
	double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
	double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);

	double lambda = L, lambdaP;
	int iterLimit = 100;
	double sigma, cosSqAlpha, sinSigma, cosSigma, cos2SigmaM;
	do {
	    double sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
	    sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
				 (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
	    if (sinSigma==0) return 0;  // co-incident points
	    cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
	    sigma = Math.atan2(sinSigma, cosSigma);
	    double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
	    cosSqAlpha = 1 - sinAlpha*sinAlpha;
	    cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
	    if (Double.isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (¤6)
	    double 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 && --iterLimit>0);

	if (iterLimit==0) return Double.NaN;  // formula failed to converge

	double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
	double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
	double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
	double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
							B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
	double s = b*A*(sigma-deltaSigma);

	return s;
    }

    /** Return distance between two points in meters, based on haversine formula */
    public static double haversineDistance(double lat1Deg, double lng1Deg, double lat2Deg, double lng2Deg) {
	double earthRadius = EARTH_RADIUS_METERS;
	double dLat = Math.toRadians(lat2Deg-lat1Deg);
	double dLng = Math.toRadians(lng2Deg-lng1Deg);
	double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
	    Math.cos(Math.toRadians(lat1Deg)) * Math.cos(Math.toRadians(lat2Deg)) *
	    Math.sin(dLng/2) * Math.sin(dLng/2);
	double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
	return earthRadius * c;
    }


    /** Return bearing (radians) from lat1,lon2 to lat2,lon2 */
    public static double bearing(double lat1Deg, double lon1Deg, double lat2Deg, double lon2Deg) {
	double f = 1/298.257223563;  // WGS-84 ellipsoid params
	double U1 = Math.atan((1-f) * Math.tan(Math.toRadians(lat1Deg)));
	double U2 = Math.atan((1-f) * Math.tan(Math.toRadians(lat2Deg)));
	double sinU1 = Math.sin(U1);
	double cosU1 = Math.cos(U1);
	double sinU2 = Math.sin(U2);
	double cosU2 = Math.cos(U2);
	double lambda = Math.toRadians(lon2Deg-lon1Deg);
	double sinLambda = Math.sin(lambda);
	double cosLambda = Math.cos(lambda);

	return Math.atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
    }



    /** Latitude and longitude (degrees) at specified epoch seconds */
    static class Location {
	long _epochSec;
	double _latitude;
	double _longitude;

	Location(long epochSec, double lat, double lon) {
	    _epochSec = epochSec;
	    _latitude = lat;
	    _longitude = lon;
	}
    }

}