function Sysgeo() { }

Sysgeo.d2r = Math.PI / 180.0;
Sysgeo.R2D = 180.0 / Math.PI;
Sysgeo.sqr = function(x) { return x * x; }

function Ellipsoid(v1, v2, code) {
    switch (code) {
    case 0: // AF
        this.a = v1;
        var f1 = v2;
        this.b = this.a * (1 - 1 / f1);
        this.e2 = (2 * f1 - 1) / (f1 * f1);
        this.e = Math.sqrt(this.e2);
        break;
    case 1: // AB
        this.a = v1;
        this.b = v2;
        this.e2 = 1 - Sysgeo.sqr(this.b / this.a); // mieux pour les spheres A = B
        this.e = Math.sqrt(this.e2);
//      var f1 = this.a / (this.a - this.b);
        break;
    case 2: // AE
        this.a = v1;
        this.e = v2;
        this.e2 = this.e * this.e;
        this.b = this.a * Math.sqrt(1 - this.e2);
        // var f1 = this.a / (this.a - this.b);
        break;
    }
    this.ep2 = Sysgeo.sqr(this.a / this.b) - 1;
}

Ellipsoid.prototype.geoTopo = function(geo) {
    var alt = 0;
    var m = geo.lon * Sysgeo.d2r;
    var lat_r = geo.lat * Sysgeo.d2r, sl = Math.sin(lat_r), cl = Math.cos(lat_r);
    var gn = this.a / Math.sqrt(1 - this.e2 * sl * sl);
    var k = (gn + alt) * cl;
    var x = k * Math.cos(m);
    var y = k * Math.sin(m);
    var z = (gn * (1 - this.e2) + alt) * sl;
    return { x: x, y: y, z: z };
};

Ellipsoid.prototype.topoGeo = function(topo) {
    var mer = 0;
    var r2 = topo.x * topo.x + topo.y * topo.y;
    if (r2 == 0.0) { // Axe Z
        var lon = 0; // convention
        var lat = topo.z < 0 ? -90 : 90;
        var alt = topo.z - this.b;
        return { lat: lat, lon: lon, alt: alt };
    }
    var lon = Math.atan2(topo.y, topo.x) * Sysgeo.R2D - mer;

    var ae2 = this.a * this.e2;
    var c1e2 = 1 - this.e2;
    var z2 = topo.z * topo.z;
    var r = Math.sqrt(r2);
    var dr = r * (1 - ae2 / Math.sqrt(r2 + z2 / c1e2));
    var u = Math.sqrt(dr * dr + z2);
    var s = topo.z / u;
    var c = dr / u;
    var s2 = s * s;
    var w2 = 1 - this.e2 * s2;
    u = ae2 / Math.sqrt(w2) / w2;
    dr = r - u * c * c * c;
    var dz = topo.z + u * c1e2 * s * s2;
    var lat = Math.atan2(dz, dr) * Sysgeo.R2D;
    u = Math.sqrt(dr * dr + dz * dz);
    s = dz / u;
    var alt = r * dr / u + topo.z * s - this.a * Math.sqrt(1 - this.e2 * s * s);
    return { lat: lat, lon: lon, alt: alt };
};

var wgs84_e     = new Ellipsoid(6378137.0, 298.257223563, 0);
var clark80_ign = new Ellipsoid(6378249.2, 6356515.0, 1);
var hayford     = new Ellipsoid(6378388.0 ,297.0, 0);

// meridiens
var paris = 2.337229167;

// la grille
function getTranslation(lat, lon) {
    var g = gr3df97a;
    if (lon < g.lon1 || lon >= g.lon2 || lat < g.lat1 || lat >= g.lat2)
        return { dx: g.dx, dy: g.dy, dz: g.dz };
    var x1 = (lon - g.lon1) / g.dlon, y1 = (lat - g.lat1) / g.dlat;
    var x0 = Math.floor(x1), y0 = Math.floor(y1);
    x1 -= x0; y1 -= y0;
    var x2 = 1 - x1, y2 = 1 - y1;
    var k1 = x2 * y2, k2 = x1 * y2, k3 = x2 * y1, k4 = x1 * y1;
    var i1 = (x0 * g.ny + y0) * 3, i2 = i1 + g.ny * 3, i3 = i1 + 3, i4 = i2 + 3;
    var dx = g.dx + (g.val[i1 + 0] * k1 + g.val[i2 + 0] * k2 + g.val[i3 + 0] * k3 + g.val[i4 + 0] * k4) * 0.001;
    var dy = g.dy + (g.val[i1 + 1] * k1 + g.val[i2 + 1] * k2 + g.val[i3 + 1] * k3 + g.val[i4 + 1] * k4) * 0.001;
    var dz = g.dz + (g.val[i1 + 2] * k1 + g.val[i2 + 2] * k2 + g.val[i3 + 2] * k3 + g.val[i4 + 2] * k4) * 0.001;
//            return "<br/>x1:" + x1 + "<br/>y1:" + y1 + "<br/>x0:" + x0 + "<br/>y0:" + y0 + "<br/>k1:" + k1 + "<br/>k2:" + k2 + "<br/>k3:" + k3 + "<br/>k4:" + k4 +
//             "<br/>kt:" + (k1 + k2 + k3 + k4) +
//             "<br/>dx:" + dx +
//             "<br/>dy:" + dy +
//             "<br/>dz:" + dz;
    return { dx: dx, dy: dy, dz: dz };
   }

function Datum(code, elli, meri, dx, dy, dz) {
    this.code = code;
    this.elli = elli;
    this.meri = meri;
    this.dx = dx;
    this.dy = dy;
    this.dz = dz;
}
wgs84 = new Datum(0, wgs84_e, 0, 0, 0, 0);
ntf = new Datum(1, clark80_ign, paris, -168, -60, +320);
var belg72 = new Datum(0, hayford,  0, -127,  79, -108);

Datum.prototype.toGeo = function(geo, datum) {
    if (this.code == 0 && datum.code == 0)
        return geo;
    if (this.code != 0) {
        alert("todo toGeo from ntf");
        return geo;
    }
    var lat = geo.lat;
    var lon = geo.lon + this.meri;
    var alt = 0;
    var topo = this.elli.geoTopo({ lat: lat, lon: lon, alt: alt });
    var trs = getTranslation(lat, lon);
    topo.x -= trs.dx;
    topo.y -= trs.dy;
    topo.z -= trs.dz;
//    alert("topo2:" + topo.x + " " + topo.y + " " + topo.z);
    geo = datum.elli.topoGeo(topo);
    geo.lon -= datum.meri;
    return geo;
}

function latIso(lat, e) {
    var esf = e * Math.sin(lat);
    return Math.log(Math.tan(Math.PI * 0.25 + lat * 0.5)) - e * 0.5 * Math.log((1 + esf) / (1 - esf));
}

function Lambert(d, m, n, c, x, y) {
    this.datum = d;
    this.m = m;
    this.n = n;
    this.c = c;
    this.x = x;
    this.y = y;
}

function LambertTangent(datum, k0, lat, lon, x0, y0) {
    var l0 = lat * Sysgeo.d2r;
    var m0 = lon * Sysgeo.d2r;
    var elli = datum.elli;
    var e = elli.e;
    var n = Math.sin(l0);
    var esf = e * n;
    var lc = latIso(l0, e);
    var k0r0 = k0 * elli.a / Math.sqrt(1 - esf * esf) * Math.cos(l0) / n;
    var c = k0r0 * Math.exp(n * lc);
    return new Lambert(datum, m0, n, c, x0, k0r0 + y0);
}

function LambertSecant(datum, lat1, lat2, lat, lon, x0, y0) {
    var m0 = lon * Sysgeo.d2r;
    var l0 = lat * Sysgeo.d2r;
    var l1 = lat1 * Sysgeo.d2r;
    var l2 = lat2 * Sysgeo.d2r;
    var elli = datum.elli;

    var e = elli.e;

    var esf1 = e * Math.sin(l1);
    var lc1 = latIso(l1, e);
    var m1 = Math.cos(l1) / Math.sqrt(1 - esf1 * esf1);

    var esf2 = e * Math.sin(l2);
    var lc2 = latIso(l2, e);
    var m2 = Math.cos(l2) / Math.sqrt(1 - esf2 * esf2);

    var n = Math.log(m1 / m2) / (lc2 - lc1);
    var c = m1 * Math.exp(n * lc1) / n * elli.a;
    var esf0 = e * Math.sin(l0);
    var lc0 = latIso(l0, e);
    var y = c * Math.exp(-n * lc0) + y0;
    var x = x0;
    return new Lambert(datum, m0, n, c, x, y);
}

Lambert.prototype.geoCarto = function(geo) {
    var e = this.datum.elli.e;
    var lc = latIso(geo.lat * Sysgeo.d2r, e);
    var r = this.c * Math.exp(-this.n * lc);
    var a = (geo.lon * Sysgeo.d2r - this.m) * this.n;
    var x = this.x + r * Math.sin(a);
    var y = this.y - r * Math.cos(a);
    return { x: x, y: y };
};

var lamb1 = LambertTangent(ntf, 0.99987734, 49.5, 0, 600000, 1200000);
var lamb3 = LambertTangent(ntf, 0.99987750, 44.1, 0, 600000, 3200000);
var lamb2e = LambertTangent(ntf, 0.99987742, 46.8, 0, 600000, 2200000);
var lamb93 = LambertSecant(wgs84, 44, 49, 46.5, 3, 700000, 6600000);
var cc42 = LambertSecant(wgs84, 41.25, 42.75, 42, 3, 1700000, 1200000);
var cc43 = LambertSecant(wgs84, 42.25, 43.75, 43, 3, 1700000, 2200000);
var cc44 = LambertSecant(wgs84, 43.25, 44.75, 44, 3, 1700000, 3200000);
var cc45 = LambertSecant(wgs84, 44.25, 45.75, 45, 3, 1700000, 4200000);
var cc46 = LambertSecant(wgs84, 45.25, 46.75, 46, 3, 1700000, 5200000);
var cc47 = LambertSecant(wgs84, 46.25, 47.75, 47, 3, 1700000, 6200000);
var cc48 = LambertSecant(wgs84, 47.25, 48.75, 48, 3, 1700000, 7200000);
var cc49 = LambertSecant(wgs84, 48.25, 49.75, 49, 3, 1700000, 8200000);
var cc50 = LambertSecant(wgs84, 49.25, 50.75, 50, 3, 1700000, 9200000);

var lamb_belge = LambertSecant(belg72, 49.83333333333333333, 51.1666666666666, 90, 4.3674866666666666, 150000.013, 5400088.438);