iseaplane.js

"use strict";
var sin = Math.sin,
    cos = Math.cos,
    tan = Math.tan,
    acos = Math.acos,
    asin = Math.asin,
    atan2 = Math.atan2,
    PI = Math.PI;

var ISEA_PLANE = Symbol('ISEA_PLANE'),
    SNYDER_POLY_ICOSAHEDRON = 6;

var constants = [
    {g: 23.80018260, G: 62.15458023, theta: 60.0, ea_w: 3.75, ea_a: 1.033, ea_b: 0.968, g_w: 5.09, g_a: 1.195, g_b: 1.0},
    {g: 20.07675127, G: 55.69063953, theta: 54.0, ea_w: 2.65, ea_a: 1.030, ea_b: 0.983, g_w: 3.59, g_a: 1.141, g_b: 1.027},
    {g: 0.0, G: 0.0, theta: 0.0, ea_w: 0.0, ea_a: 0.0, ea_b: 0.0, g_w: 0.0, g_a: 0.0, g_b: 0.0},
    {g: 0.0, G: 0.0, theta: 0.0, ea_w: 0.0, ea_a: 0.0, ea_b: 0.0, g_w: 0.0, g_a: 0.0, g_b: 0.0},
    {g: 0.0, G: 0.0, theta: 0.0, ea_w: 0.0, ea_a: 0.0, ea_b: 0.0, g_w: 0.0, g_a: 0.0, g_b: 0.0},
    {g: 0.0, G: 0.0, theta: 0.0, ea_w: 0.0, ea_a: 0.0, ea_b: 0.0, g_w: 0.0, g_a: 0.0, g_b: 0.0},
    {g: 37.37736814, G: 36.0, theta: 30.0, ea_w: 17.27, ea_a: 1.163, ea_b: 0.860, g_w: 13.14, g_a: 1.584, g_b: 1.0}
];


var TABLE_G = 0.6615845383,
    TABLE_H = 0.1909830056,
    RPRIME = 0.91038328153090290025,
    ISEA_STD_LAT = 1.01722196792335072101,
    ISEA_STD_LON = 0.19634954084936207740,
    PRECISION = 0.0000000000005;

var E = 52.62263186,
    F = 10.81231696,
    DEG60 = 1.04719755119659774614,
    DEG120 = 2.09439510239319549229,
    DEG72 = 1.25663706143591729537,
    DEG90 = 1.57079632679489661922,
    DEG144 = 2.51327412287183459075,
    DEG36 = 0.62831853071795864768,
    DEG108 = 1.88495559215387594306,
    DEG180 = PI,
    ISEA_SCALE = 0.8301572857837594396028083,
    V_LAT = 0.46364760899944494524,
    RAD2DEG = (180 / PI),
    DEG2RAD = (PI / 180),
    E_RAD = 0.91843818702186776133,
    F_RAD = 0.18871053072122403508;

var tri_v1 = [0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11];

var vertex = [
    [0.0, DEG90],
    [DEG180, V_LAT],
    [-DEG108, V_LAT],
    [-DEG36, V_LAT],
    [DEG36, V_LAT],
    [DEG108, V_LAT],
    [-DEG144, -V_LAT],
    [-DEG72, -V_LAT],
    [0.0, -V_LAT],
    [DEG72, -V_LAT],
    [DEG144, -V_LAT],
    [0.0, -DEG90]
];

var icostriangles = [
    [0.0, 0.0],
    [-DEG144, E_RAD],
    [-DEG72, E_RAD],
    [0.0, E_RAD],
    [DEG72, E_RAD],
    [DEG144, E_RAD],
    [-DEG144, F_RAD],
    [-DEG72, F_RAD],
    [0.0, F_RAD],
    [DEG72, F_RAD],
    [DEG144, F_RAD],
    [-DEG108, -F_RAD],
    [-DEG36, -F_RAD],
    [DEG36, -F_RAD],
    [DEG108, -F_RAD],
    [DEG180, -F_RAD],
    [-DEG108, -E_RAD],
    [-DEG36, -E_RAD],
    [DEG36, -E_RAD],
    [DEG108, -E_RAD],
    [DEG180, -E_RAD]
];

function isea_triangle_xy(triangle){
    triangle = (triangle - 1) % 20;
    var c = {};
    c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
    if (triangle > 9) {
        c.x += TABLE_G;
    }
    let t = Math.floor(triangle / 5);
    if(t == 0){
        c.y = 5.0 * TABLE_H;
    } else if (t == 1){
        c.y = TABLE_H;
    } else if (t == 2){
        c.y = -TABLE_H;
    } else if (t == 3){
        c.y = -5.0 * TABLE_H;
    } else {
        c.y = -5.0 * TABLE_H;
        console.log("Error; t = ", t);
        return;
    }
    c.x *= RPRIME;
    c.y *= RPRIME;
    return c;
}

function isea_rotate(pt, degrees){
    var x, y,
        rad = -degrees * PI / 180.0;

    while (rad >= 2.0 * PI) rad -= 2.0 * PI;
    while (rad <= -2.0 * PI) rad += 2.0 * PI;

    x = pt.x * cos(rad) + pt.y * sin(rad);
    y = -pt.x * sin(rad) + pt.y * cos(rad);

    pt.x = x;
    pt.y = y;
}

function _downtri(tri){ return (((tri - 1) / 5) % 2 == 1); }

function isea_tri_plane(tri, pt, radius){
    var tc = {};
    if(_downtri(tri)){
        isea_rotate(pt, 180);
    }
    tc = isea_triangle_xy(tri);
    tc.x *= radius;
    tc.y *= radius;
    pt.x += tc.x;
    pt.y += tc.y;
    return tri;
}

function az_adjustment(triangle){
    var v = vertex[tri_v1[triangle]],
        c = icostriangles[triangle];

    v = {lon: v[0], lat: v[1]};
    c = {lon: c[0], lat: c[1]};
    return atan2(
        cos(v.lat) * sin(v.lon - c.lon),
        cos(c.lat) * sin(v.lat) - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon)
      );
}

function sph_azimuth(f_lon, f_lat, t_lon, t_lat){
    return atan2(
        cos(t_lat) * sin(t_lon - f_lon),
        cos(f_lat) * sin(t_lat) - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
      );
}

function isea_snyder_forward(ll, out){
    var i, g, G, theta, q, H, Ag, Azprime, Az, dprime, f, rho, x, y;
    var cot_theta, tan_g, az_offset, Az_adjust_multiples, c;
    c = constants[SNYDER_POLY_ICOSAHEDRON];
    theta = c.theta * DEG2RAD;
    g = c.g * DEG2RAD;
    G = c.G * DEG2RAD;
    for(i=1; i <= 20; i++){
        var center = icostriangles[i], z;
        center = {lon: center[0], lat: center[1]};
        z = acos(sin(center.lat) * sin(ll.lat) + cos(center.lat) * cos(ll.lat) * cos(ll.lon - center.lon));
        if (z > g + 0.000005){
            continue;
        }
        Az = sph_azimuth(center.lon, center.lat, ll.lon, ll.lat);
        az_offset = az_adjustment(i);
        Az -= az_offset;
        if (Az < 0.0) {
            Az += 2.0 * PI;
        }
        Az_adjust_multiples = 0;
        while (Az < 0.0) {
            Az += DEG120;
            Az_adjust_multiples--;
        }
        while (Az > DEG120 + 0.000005) {
            Az -= DEG120;
            Az_adjust_multiples++;
        }
        cot_theta = 1.0 / tan(theta);
        tan_g = tan(g);
        q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);

        if (z > q + 0.000005) {
            continue;
        }
        H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
        Ag = Az + G + H - DEG180;
        Azprime = atan2(2.0 * Ag, RPRIME * RPRIME * tan_g * tan_g - 2.0 * Ag * cot_theta);
        dprime = RPRIME * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
        f = dprime / (2.0 * RPRIME * sin(q / 2.0));
        rho = 2.0 * RPRIME * f * sin(z / 2.0);
        Azprime += DEG120 * Az_adjust_multiples;
        x = rho * sin(Azprime);
        y = rho * cos(Azprime);
        out.x = x;
        out.y = y;
        return i;
    }
    console.log("Error - Not supposed to happend");
}

function snyder_ctran(np, pt){
    var npt = {},
        alpha, phi, lambda, lambda0, beta,
        lambdap, phip, sin_phip, lp_b, cos_p, sin_a;

    phi = pt.lat; lambda = pt.lon;
    alpha = np.lat; beta = np.lon;
    lambda0 = beta;
    cos_p = cos(phi);
    sin_a = sin(alpha);
    sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);

    lp_b = atan2(cos_p * sin(lambda - lambda0),
       (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
    lp_b = atan2(cos_p * sin(lambda - lambda0),
       (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));

    lambdap = lp_b + beta;
    lp_b = atan2(cos_p * sin(lambda - lambda0),
       (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));

    lambdap = lp_b + beta;

    lambdap = lp_b + beta;

    lambdap = lambdap % (2 * PI);
    while (lambdap > PI)
        lambdap -= 2 * PI;
    while (lambdap < -PI)
        lambdap += 2 * PI;

    phip = asin(sin_phip);

    npt.lat = phip;
    npt.lon = lambdap;

    return npt;
}

function isea_ctran(np, pt, lon0){
    var npt = {};
    np.lon += PI;
    npt = snyder_ctran(np, pt);
    np.lon -= PI;
    npt.lon -= (PI - lon0 + np.lon);
    npt.lon += PI;
    npt.lon = npt.lon %  (2 * PI);
    while (npt.lon > PI)
        npt.lon -= 2 * PI;
    while (npt.lon < -PI)
        npt.lon += 2 * PI;
    return npt;
}

function isea_grid_init(g){
    if(!g) return 0;
    g['polyhedron'] = 20;
    g['o_lat'] = ISEA_STD_LAT;
    g['o_lon'] = ISEA_STD_LON;
    g['o_az'] = 0;
    g['aperture'] = 4;
    g['resolution'] = 6;
    g['radius'] = 1.0;
    g['topology'] = 6;
    return 1;
}

function isea_orient_isea(g){
    if(!g) return 0;
    g['o_lat'] = ISEA_STD_LAT;
    g['o_lon'] = ISEA_STD_LON;
    g['o_az'] = 0;
    return 1;
}

function isea_orient_pole(g){
    if(!g) return 0;
    g['o_lat'] = PI / 2;
    g['o_lon'] = 0;
    g['o_az'] = 0;
    return 1;
}

function isea_transform(g, _in, _out){
    var i, pole = {}, tri;
    pole.lat = g.o_lat;
    pole.lon = g.o_lon;
    i = isea_ctran(pole, _in, g['o_az']);
    tri = isea_snyder_forward(i, _out);
    _out.x *= g.radius;
    _out.y *= g.radius;
    g.triangle = tri;
    return tri;
}

var isea_forward = function(params, _in){
    var tri, _out = {}, coord;
    tri = isea_transform(params, _in, _out);
    if(params['output'] == ISEA_PLANE){
        isea_tri_plane(tri, _out, params['radius']);
        return _out;
    }
}

function IcosahedralSnyderEqualArea(orient="isea"){
    var params = {};
    isea_grid_init(params);
    params['output'] = ISEA_PLANE;
    // params['resolution'] = 4;
    // params['aperture'] = 3;
    if(orient == "isea"){
        isea_orient_isea(params);
    } else if (orient == "pole"){
        isea_orient_pole(params);
    }

    return function(lambda, phi){
        var out = isea_forward(params, {lon: lambda, lat: phi});
        return [out.x, out.y];
    };
}