Lambert projection for Modest Maps

lambert.js

var LambertProjection = function(lat0, lon0, lat1, lat2, xmin, ymin, xmax, ymax)
{
    var pi = Math.PI, ln = Math.log, pow = Math.pow,
        sin = Math.sin, cos = Math.cos, tan = Math.tan,
        atan = Math.atan, sqrt = Math.sqrt;
    
    function sec(t) { return 1 / cos(t); }
    function cot(t) { return 1 / tan(t); }

    function deg2rad(deg) { return pi * deg / 180; }

    var phi0 = deg2rad(lat0),
        lam0 = deg2rad(lon0),
        phi1 = deg2rad(lat1),
        phi2 = deg2rad(lat2);
    
    this.lam0 = lam0;
    
    // Equation 4
    this.n = ln(cos(phi1) * sec(phi2));
    this.n /= ln(tan(pi/4 + phi2/2) * cot(pi/4 + phi1/2));

    // Equation 3
    this.F = cos(phi1) * pow(tan(pi/4 + phi1/2), this.n) / this.n;
    
    // Equation 6
    this.P0 = this.F * pow(cot(pi/4 + phi0/2), this.n);
    
    //

    var tx = 1 / (xmax - xmin),
        dx = .5 - (xmin/2 + xmax/2) * tx;
    
    var ty = -tx,
        dy = .5 - (ymin/2 + ymax/2) * ty;
    
    var t = new MM.Transformation(tx, 0, dx, 0, ty, dy);
    
    MM.Projection.call(this, 0, t);
};

LambertProjection.prototype = {
    rawProject: function(point)
    {
        var pi = Math.PI, ln = Math.log, pow = Math.pow,
            sin = Math.sin, cos = Math.cos, tan = Math.tan,
            atan = Math.atan, sqrt = Math.sqrt;
        
        function cot(t) { return 1 / tan(t); }

        var lam = point.x, // longitude already in radians
            phi = point.y; // latitude already in radians
        
        // Equation 5
        var P = this.F * pow(cot(pi/4 + phi/2), this.n);
        
        // Equations 1 & 2
        var x = P * sin(this.n * (lam - this.lam0)),
            y = this.P0 - P * cos(this.n * (lam - this.lam0));
        
        return {x: x * 6378137., y: y * 6378137.};
    },

    rawUnproject: function(point)
    {
        var pi = Math.PI, ln = Math.log, pow = Math.pow,
            sin = Math.sin, cos = Math.cos, tan = Math.tan,
            atan = Math.atan, sqrt = Math.sqrt;
        
        function sgn(x)
        {
            if(x >= 0) {
                return 1;

            } else{
                return -1;
            }
        }
        
        var x = point.x / 6378137.,
            y = point.y / 6378137.;
        
        // Equation 9
        var P = sgn(this.n) * sqrt(pow(x, 2) + pow((this.P0 - y), 2));
        
        // Equation 10
        var theta = atan(x / (this.P0 - y));
        
        // Equation 7
        var phi = 2 * atan(pow(this.F/P, 1/this.n)) - pi/2;
        
        // Equation 8
        var lam = this.lam0 + theta/this.n;
        
        return {x: lam, y: phi};
    }
};

MM.extend(LambertProjection, MM.Projection);

use-it.js

// adjust for latitude.
var scale = Math.cos(Math.PI * 37.5/180);

// echo the usual +/- 20m extent of spherical mercator.
var xmin = scale * -6378137 * Math.PI,
    ymin = scale * -6378137 * Math.PI,
    xmax = scale * 6378137 * Math.PI,
    ymax = scale * 6378137 * Math.PI;

var lcc = new LambertProjection(37.5, -96, 29.5, 45.5, xmin, ymin, xmax, ymax);

var loc = {y: Math.PI * 48.38544/180, x: Math.PI * -124.72916/180};
var proj = lcc.rawProject(loc);
var loc2 = lcc.rawUnproject(proj);

console.log([loc, proj, loc2]);

var loc3 = {lat: 48.38544, lon: -124.72916};
var coord = lcc.locationCoordinate(loc3);
var loc4 = lcc.coordinateLocation(coord);

console.log([loc3, coord, loc4]);