/********************************************************************
 * Map related functions
 *
 * Example:
    import map, display;
    G = map(".\\data\\gshhs_c.b", 110, 180, 10, 10, 80, 10);
    show(G.render);
 *
 ********************************************************************/
 
DEG_TO_RAD = 0.017453292519943;     // degree to radius


/********************************************************************
 * Unprojected Map
 *
 * Input:
 *      gshhs -- GSHHS data file
 *      lon1, lon2, dlon  -- longitude range and tickmark step.
 *      lat1, lat2, dlat  -- latitude range and tickmark step.
 *
 * Output:
 *      An ctalk array that contains render, scene, node, plot, coast,
 *      grid, font, xaxis, yaxis, and frame.
 *
 ********************************************************************/

load("zegraph.dll", "matrix.dll");


__scale = 1;

function map(gshhs, lon1=0, lon2=360, dlon=30, lat1=-90, lat2=90, dlat=30, scale=1)
{
    assert(lon1 < lon2, dlon > 0, lat1 < lat2, dlat > 0);
    
    global:__scale = scale;
    w = scale * 600.0;
    h = w * (lat2 - lat1) / (lon2 - lon1);

    coast = coast(gshhs, lon1, lon2, lat1, lat2, 0);
    grid = grid(lon1, lon2, dlon, lat1, lat2, dlat, 0);
    
    G = zegraph("render", "scene", "node", "plot", "font");
    G.coast = coast;
    G.grid = grid;
    
    render = G.render;
    render.size(w, h);
    scene = G.scene;
    scene.viewport(0, 0, w, h);
    node = G.node;
    plot = G.plot;    
    render.add(scene);
    scene.root(node);
    node.add(plot);
    plot.add(coast, grid);
    plot.font(G.font, 10*scale);
    
    axis = plot.xaxis(0, -1, 0);
    G.xaxis = axis;
    axis.range(lon1, lon2);
    axis.tickmarks(lon1, dlon, 0);
    axis.tickdigits(0, false);
    axis.title("Longitude (deg-E)");

    axis = plot.yaxis(-1, 0, 0);
    G.yaxis = axis;
    axis.range(lat1, lat2);
    axis.tickmarks(lat1, dlat, 0);
    axis.tickdigits(0, false);
    axis.title("Latitude (deg-N)");
          
    G.frame = zegraph("line");
    frame = G.frame;
    xy = zegraph("vertex");
    frame.solid(1.5*scale);
    frame.vertex(xy);
    frame.type("loop");
    xy.add(lon1, lat1, 0, lon2, lat1, 0, lon2, lat2, 0, lon1, lat2, 0);
    plot.add(frame);

    return G;
}


/********************************************************************
 * Extract GSHHS data for 2D or 3D plot
 ********************************************************************/
function coast(gshhs, lon1=0, lon2=360, lat1=-90, lat2=90, z=0)
{
    assert(lon1 < lon2, lat1 < lat2);
    
    x1 = lon1;
    x2 = lon2;
    if (lon1 < 0) {
        x1 += 180;
        x2 += 180;
    }
    
    data = matrix("double");
    data.readgshhs(gshhs, "land", x1, x2, lat1, lat2);
    data.insert(2, z);
    
    if (lon1 < 0) {
        a = data[null,0];
        a -= 180;
        data[null,0] = a;
    }

    xyz = zegraph("vertex");
    p = data.ptr();
    xyz.add(p[0], p[1]);
        
    coast = zegraph("line");
    coast.vertex(xyz);
    coast.type("lines");
    
    return coast;
}


/********************************************************************
 * Make map grid
 ********************************************************************/
function grid(lon1=0, lon2=360, dlon=30, lat1=-90, lat2=90, dlat=30, z=0)
{
    assert(lon1 < lon2, dlon > 0, lat1 < lat2, dlat > 0);
  
    grid = zegraph("line");
    xyz = zegraph("vertex");
    grid.dot(1);
    grid.vertex(xyz);
    grid.type("lines");
    grid.color(0.5, 0.5, 0.5);
    
    for (a = lon1 + dlon;  a < lon2; a += dlon) {
        xyz.add(a, lat1, z,  a, lat2, z);
    }
    for (a = lat1 + dlat;  a < lat2; a += dlat) {
        xyz.add(lon1, a, z,  lon2, a, z);
    }
    
    return grid;
}


/********************************************************************
 * Sphere projection
 ********************************************************************/

function spheric(lon, lat)
{
    lat = DEG_TO_RAD * lat;
    lon = DEG_TO_RAD * lon;
    z = sin(lat);
    x = cos(lat) * cos(lon);
    y = cos(lat) * sin(lon);
    return [x, y, z];
}

/********************************************************************
 * Aitoff projection.
 *
 * Based on the equations
 *      x = 2 * a * cos(lat) * sin(0.5 *lon) / sin(a)
 *      y = a * sin(lat) / sin(a)
 *      cos(a) = cos(lat) * cos(0.5 * lon)
 *       
 * of Cartographic Projection Procedures
 * by G.I. Evenden, September 24, 1995
 *  
 * lon and lat are expected to vary from -180 to 180 and
 * from -90 to 90, respectively. Resulted lon varies from
 * -3.06 to 3.10 and lat from -1.48 to 1.54.
 ********************************************************************/

function aitoff(lon, lat)
{
    lon = DEG_TO_RAD * lon;
    lon /= 2.0;
    lat = DEG_TO_RAD * lat;
    a = acos(cos(lat) * cos(lon));
    a /= (1.e-5 + sin(a));
    x = 2.0 * a * cos(lat) * sin(lon);
    y = a * sin(lat);
    return [x, y];
}

/********************************************************************
 * Denoyer projection.
 *    
 * Based on the equations:
 *      x = lon * cos((0.95 - lon/12 + lon^3/600) * lat)
 *       y = lat
 * of Cartographic Projection Procedures by
 * G.I. Evenden, September 24, 1995
 *  
 * lon and lat are expected to vary from -180 to 180 and
 * from -90 to 90, respectively. Resulted lon varies from
 * -3.07 to 3.10 and lat from -1.37 to 1.46.
 ********************************************************************/

function denoyer(lon, lat)
{
    lon = DEG_TO_RAD * lon;
    lat = DEG_TO_RAD * lat;
    lon *= cos((0.95 - abs(lon) / 12.0 + cos(pow(lon, 3)) / 600.0) * lat);  
    return [lon, lat];
}

/********************************************************************
 * Eckert V projection.
 *
 * Based on the equations
 *      x = lon * (1 + cos(lat)) / sqrt(2 + pi)
 *      y = 2 * lat / sqrt(2 + pi)
 * of Cartographic Projection Procedures
 * by G.I. Evenden, September 24, 1995
 *   
 * lon and lat are expected to vary from -180 to 180 and
 * from -90 to 90, respectively. Resulted lon varies from
 * -2.71 to 2.74 and lat from -1.21 to 1.28.
 ********************************************************************/

function echert(lon, lat)
{
    lon = DEG_TO_RAD * lon;
    lat = DEG_TO_RAD * lat;
    lon *= (1 + cos(lat)) *  0.441012772;
    lat *= 0.882025544;
    return [lon, lat];
}

/********************************************************************
 * Wagner II projection.
 *
 * Based on the equations:
 *      x = 0.92483 * lon * cos(theta)
 *      y = 1.38725 * theta
 *      sin(theta) = 0.88022 * sin(0.8855 * lat)
 * of Cartographic Projection Procedures
 * by G.I. Evenden, September 24, 1995
 *   
 * lon and lat are expected to vary from -180 to 180 and
 * from -90 to 90, respectively. Resulted lon varies from
 * -2.84 to 2.84 and lat from -1.34 to 1.34.
 ********************************************************************************/

function wagner(lon, lat)
{
    lon = DEG_TO_RAD * lon;
    lat = DEG_TO_RAD * lat;
    theta = asin(0.88022 * sin(0.8855 * lat));
    lon *= 0.92483 * cos(theta);
    lat = 1.38725 * theta;
    return [lon, lat];
}