proj4.js

// Transliteration of a couple functions from PROJ.4 to JavaScript.
// I have absolutely no idea how any of this works.

/** Some constants that PROJ uses for the WGS84 ellipse projection.*/
var G_A = 6378137;
var G_ONEF = 0.99664718933525254;
var G_FLAT4 = 0.00083820266618686579;
var G_FLAT64 = 1.756459274006926e-07;
var METERS_TO_MILES = 0.00062136994949494966;
var DTOL = 1e-12;
var SPI = 3.14159265359;


// Given two latitude/longitude tuples, returns the distance between them using the
// WGS84 projection (which I believe is the IATA standard).
function distanceInMiles(lat1, lon1, lat2, lon2) {
     // This does...something.  Apparently, it's important.
    function adjlon(lon) {
    	if(Math.abs(lon) <= SPI) return lon;
    	lon += SPI;  /* adjust to 0..2pi rad */
    	lon -= 2 * SPI * Math.floor(lon / (2 * SPI)); /* remove integral # of 'revolutions'*/
    	lon -= SPI;  /* adjust back to -pi..pi rad */
    	return lon;
    }
    
    // Convert to radians
	lat1 = lat1 * (SPI / 180);
	lat2 = lat2 * (SPI / 180);
	lon1 = lon1 * (SPI / 180);
	lon2 = lon2 * (SPI / 180);
	var th1 = Math.atan(G_ONEF * Math.tan(lat1));
	var th2 = Math.atan(G_ONEF * Math.tan(lat2));
	var thm = .5 * (th1 + th2);
	var dthm = .5 * (th2 - th1);
	var dlam = 0;
	var dlamm = .5 * ( dlam = adjlon(lon2 - lon1) );
    var geod_S = 0;
	if(Math.abs(dlam) < DTOL && Math.abs(dthm) < DTOL) {
		geod_S = 0.;
		return NaN;
	}
	var sindlamm = Math.sin(dlamm);
	var costhm = Math.cos(thm);
	var sinthm = Math.sin(thm);
	var cosdthm = Math.cos(dthm);
	var sindthm = Math.sin(dthm);
	var L = sindthm * sindthm + (cosdthm * cosdthm - sinthm * sinthm) * sindlamm * sindlamm;
	var cosd = 0;
	var d = Math.acos(cosd = 1 - L - L);
	// Holy crap!
	var E = cosd + cosd;
	var sind = Math.sin( d );
	var Y = sinthm * cosdthm;
	Y *= (Y + Y) / (1. - L);
	var T = sindthm * costhm;
	T *= (T + T) / L;
	var X = Y + T;
	Y -= T;
	T = d / sind;
	var D = 4. * T * T;
	var A = D * E;
	var B = D + D;
	geod_S = G_A * sind * (T - G_FLAT4 * (T * X - Y) +
        G_FLAT64 * (X * (A + (T - .5 * (A - E)) * X) -
		Y * (B + E * Y) + D * X * Y));
	return geod_S * METERS_TO_MILES;
}