toUtm33.js

/**
Simple script to transform coordinates from WGS84 Geographic to UTM Zone 33N.
Based on Holsen.js, which in turn is based on "Holsens småprogrammer"
License: MIT
Written by: Atle F. Sveen
**/

var toUtm33 = (function () {
    'use strict';

    //make these easier to type/read
    var cos = Math.cos;
    var sin = Math.sin;
    var tan = Math.tan;
    var pow = Math.pow;
    var sqrt = Math.sqrt;

    //Ellipsoid constants for WGS84-ellipsoid
    var ellipsoid = {
        a: 6378137.000,
        f: 1 / 298.257223563
    };
    ellipsoid.b = -1 * ((ellipsoid.f * ellipsoid.a) - ellipsoid.a);

    //Values specific for Zone 33N
    var coordsystem = {
        scaleFactor: 0.9996,
        falseEasting: 500000,
        falseNorthing: 0,
        centralMeridian: 15,
        latitudeOfOrigin: 0
    };

    function toRad(deg) {
        return deg / (180 / Math.PI);
    }

    function ellipsoidParams(ellipsoid) {
        var n = (ellipsoid.a - ellipsoid.b) / (ellipsoid.a + ellipsoid.b);
        return {
            n: n,
            k: ellipsoid.a / (n + 1),
            k1: 1 + n * (n / 4) + pow(n, 4) / 64,
            k2: (n - n * n * n / 8) * 3,
            k3: (n * n - pow(n, 4) / 4) * 15 / 8
        };
    }

    function getGeodeticArch(lon1, lon2, ell) {
        var db = lon1 - lon2;
        var bm = lon1 - db / 2;
        return (ell.k * (ell.k1 * db - ell.k2 * cos(2 * bm) * sin(db) +
            ell.k3 * cos(4 * bm) * sin(2 * db)) -
            ell.k * (pow(ell.n, 3) * cos(6 * bm) * sin(3 * db) * 35 / 24) +
            ell.k * (pow(ell.n, 4) * cos(8 * bm) * sin(4 * db) * 315) / 256);
    }

    return function (lon, lat) {
        lon = toRad(lon) - toRad(coordsystem.centralMeridian);
        lat = toRad(lat);

        var a = ellipsoid.a;
        var b = ellipsoid.b;
        var et = (pow(a, 2) - pow(b, 2)) / pow(b, 2) * pow(cos(lat), 2);
        var n1 = pow(a, 2) / (sqrt(1 + et) * b);
        var t = tan(lat);
        var a1 = n1 * cos(lat);
        var a2 = -(n1 * t * pow(cos(lat), 2)) / 2.0;
        var a3 = -(n1 * pow(cos(lat), 3)) * (1 - pow(t, 2) + et) / 6;
        var a4 = n1 * t * pow(cos(lat), 4) * (5 - pow(t, 2) + 9 * et + 4 * pow(et, 2)) / 24;
        var a5 = n1 * pow(cos(lat), 5) * (5 - 18 * pow(t, 2) + pow(t, 4) + 14 * et - 58 * et * pow(t, 2)) / 120;
        var a6 = n1 * t * pow(cos(lat), 6) * (61 - 58 * pow(t, 2) + pow(t, 4) + 270 * et - 330 * et * pow(t, 2)) / 720;

        var north = -a2 * pow(lon, 2) + a4 * pow(lon, 4) + a6 * pow(lon, 6);

        var east = a1 * lon - a3 * pow(lon, 3) + a5 * pow(lon, 5);
        north = north + getGeodeticArch(
            lat,
            toRad(coordsystem.latitudeOfOrigin),
            ellipsoidParams(ellipsoid)
        );

        //apply UTM "tricks"
        return {
            y: north * coordsystem.scaleFactor + coordsystem.falseNorthing,
            x: east * coordsystem.scaleFactor + coordsystem.falseEasting
        };
    };
}());