/********************************************************************
 * Functions for making globe
 *
 * Example:
    import globe, display;
    G = globe_map(200, 10, "data\\gshhs_c.b");
    animate(G.render, G.node, G.rotatex, G.rotatez);
 *
 ********************************************************************/

load("zegraph.dll", "matrix.dll");
 
DEG_TO_RAD = 0.017453292519943;     // degree to radius

__R = null;

/********************************************************************
 * Callback functions for making globe surface
 ********************************************************************/
function vertex_func(u, v)
{
    u = DEG_TO_RAD * u;     // in case of a matrix, this prevent from modifying the original copy.
    v = DEG_TO_RAD * v;
    z = __R * sin(v);
    r = __R * cos(v);
    x = r * cos(u);
    y = r * sin(u);
    return [x, y, z];
}

function normal_func(u, v)
{
    u = DEG_TO_RAD * u;
    v = DEG_TO_RAD * v;
    z = sin(v);
    r = cos(v);
    x = r * cos(u);
    y = r * sin(u);
    return [x, y, z];
}

/********************************************************************
 * grid lines
 ********************************************************************/
function globe_grid(r, grid)
{
    grid = integer(grid);
    global:__R = r;
    r = 1.001 * r;
    
	line = zegraph("line");
    xyz = zegraph("vertex");
    nor = zegraph("vertex");
    line:type("lines");
    line:color(0.4, 0.4, 0.4);
    line:vertex(xyz);
    line:normal(nor);
    grid = integer(grid);

    for (lon = 0.0; lon < 360; lon += grid) {
        lat = -90 + grid;
		a = normal_func(lon, lat);
		lat2 = 90.0 - grid;
		for (lat += 2.0; lat <= lat2; lat += 2.0) {
			b = normal_func(lon, lat);
		    xyz:add(r*a[0], r*a[1], r*a[2],  r*b[0], r*b[1], r*b[2]);
		    nor:add(a[0], a[1], a[2],  b[0], b[1], b[2]);
		    a = b;
		}
    }
    lat2 = 90 - grid;
    for (lat = -90.0 + grid; lat <= lat2; lat += grid) {
	    lon = 0;
		a = normal_func(lon, lat);
		for (lon += 2.0; lon <= 360.0; lon += 2.0) {
			b = normal_func(lon, lat);
		    xyz:add(r*a[0], r*a[1], r*a[2],  r*b[0], r*b[1], r*b[2]);
		    nor:add(a[0], a[1], a[2],  b[0], b[1], b[2]);
		    a = b;
		}
    }

    return line;
}

/********************************************************************
 * coastlines
 ********************************************************************/
function globe_coast(r, fname)
{
    r = real(r);
    global:__R = r;
    
    line = zegraph("line");
    xyz = zegraph("vertex");
    nor = zegraph("vertex");
    line:type("lines");
    line:color(1, 1, 1);
    line:vertex(xyz);
    line:normal(nor);
    
    data = matrix("double");
    data:readgshhs(fname, "land");
    lon = data[null,0];
    lat = data[null,1];

    a = normal_func(lon, lat);
    data[null,0] = a[0];
    data[null,1] = a[1];
    data:insert(2, a[2]);

    p = data:ptr();
    nor:add(p[0], p[1]);
    data *= 1.001*r;
    xyz:add(p[0], p[1]);
    
    return line;
}

/********************************************************************
 * (lon, lat)  to (x, y, z) projection. Return an array
 ********************************************************************/
function globe_projection(r, lon, lat)
{
	global:__R = r;
	return vertex_func(lon, lat);
}

/********************************************************************
 * globe with coastlines and grids
 *
 * Input:
 *		r     -- radius
 *      grid  -- grid of longitude and latitude (in degree)
 *		fname -- file name of GSHHS data
 *
 * Output:
 *      Array containing objects of render, scene, root (node), node,
 *		light, grid (line), and coast (line) for further manipulation.
 *
 ********************************************************************/
function globe_map(r, grid, fname)
{
    assert(r > 0 & grid > 0);
    
    global:__R = r;
    
    G = zegraph("render", "scene", "node",  "light");
    root = zegraph("node");
    G.root = root;
    G.rotatex = -60;
    G.rotatez = -90;
    
    G.grid = globe_grid(r, grid);
    G.coast = globe_coast(r, fname);
    
    sphere = zegraph("polygon");
    G.spere = sphere;
    sphere:sphere(r, 5);
    sphere:color(0, 0, 0.5);

    w = integer(3 * r);
    
    render = G.render;
    render:size(w, w);
    scene = G.scene;
    scene:viewport(0, 0, w, w);
    render:add(scene);
    scene:root(root);
    
    node = G.node;
    node:add(sphere, G.grid, G.coast);
    node:rotatez(G.rotatez);
    node:rotatex(G.rotatex);
    
    light = G.light;
    light:position(0.25, 0.25, 1);
    light:ambient(0.5, 0.5, 0.5);
    
    root:add(light, node);
    
    return G;
}

/********************************************************************
 * globe with texture map
 *
 * Input:
 *		r     -- radius
 *      grid  -- grid of longitude and latitude (in degree)
 *		fname -- file name of GSHHS data
 *
 * Output:
 *      Array containing objects of render, scene, root (node), node,
 *		light, grid (line), and coast (line) for further manipulation.
 *
 ********************************************************************/
function globe_texture(r, grid, fname)
{
    assert(r > 0 & grid > 0);
    
    r = real(r);
    global:__R = r;

    G = zegraph("render", "scene", "node", "light");
    G.root = zegraph("node");
    G.grid = globe_grid(r, grid);
    G.rotatex = -60;
    G.rotatez = -90;
    
    w = integer(3 * r);
    
    render = G.render;
    render:size(w, w);
    root = G.root;
    scene = G.scene;
    scene:viewport(0, 0, w, w);
    render:add(scene);
    scene:root(root);
    node = G.node;
    light = G.light;
    node:rotatez(G.rotatez);
    node:rotatex(G.rotatex);
//    node:add(G.grid);
    light:position(0.25, 0.25, 1);
    light:ambient(0.5, 0.5, 0.5);
    root:add(light, node);

    node2 =  zegraph("node");
    texture = zegraph("texture");
    polygon = zegraph("polygon");
    node2:add(texture, polygon);
    node2:closed(true);
    node:add(node2);
    
    texture:image(fname, -1, -1, -1);
    
    polygon:sphere(r, 5, true);
    
    return G;
}

/********************************************************************
 * Greate circle distance on unit globe
 *
 * Input:
 *      lon1, lat1 -- longitude/latitude of start point (in degree)
 *      lon2, lat2 -- longitude/latitude of end point (in degree)
 *
 * Output:
 *      The distance
 *
 ********************************************************************/
function globe_length(lon1, lat1, lon2, lat2)
{
	lon1 = DEG_TO_RAD * lon1;
	lat1 = DEG_TO_RAD * lat1;
	lon2 = DEG_TO_RAD * lon2;
	lat2 = DEG_TO_RAD * lat2;
	return acos(cos(lat1)*cos(lat2)*cos(lon1-lon2) + sin(lat1)*sin(lat2));
}

/********************************************************************
 * Create a line object to be drawn as an arc on a sphere surface.
 * The arc line should be less then 2 on a unit sphere.
 *
 * Input:
 *      x1, y1, z1 -- start coordinate on a unit sphere
 *      x2, y2, z2 -- end coordinate on a unit sphere
 *
 * Output:
 *      A line object. The interval between vertices is about the
 *      langth of a 1-degree arc.
 *
 ********************************************************************/
function globe_lxyz(r, x1, y1, z1, x2, y2, z2)
{
    r = real(r);
    
	line = zegraph("line");
	vertex = zegraph("vertex");
	normal = zegraph("vertex");
	line:type("strip");
	line:vertex(vertex);
	line:normal(normal);
	
	dx = x2 - x1;
	dy = y2 - y1;
	dz = z2 - z1;
	d = sqrt(dx*dx + dy*dy + dz*dz);
	n = integer(d*180.0);
	
	vertex:add(x1*r, y1*r, z1*r);
	normal:add(x1, y1, z1);
   	if (n < 1) {
		vertex:add(x2*r, y2*r, z2*r);
		normal:add(x2, y2, z2);
		return line;
	}
	
	dx /= n;
	dy /= n;
	dz /= n;
	
	for (i = 1; i < n; i++) {
		x = x1 + i*dx;
		y = y1 + i*dy;
		z = z1 + i*dz;
		d = sqrt(x*x + y*y + z*z);
		x /= d;
		y /= d;
		z /= d;
		vertex:add(x*r, y*r, z*r);
		normal:add(x, y, z);
	}

	vertex:add(x2*r, y2*r, z2*r);
	normal:add(x2, y2, z2);
	
    return line;
}

/********************************************************************
 * Create a line object to be drawn as an arc on a sphere surface.
 * The arc line should be less then 2 on a unit sphere.
 *
 * Input:
 *      lon1, lat1 -- longitude and latitude of the start point (in degree)
 *      lon2, lat2 -- longitude and latitude of the end point (in degree)
 *      r          -- radius of sphere
 *
 * Output:
 *      A line object. The interval between vertices is about the langth
 *      of a 1-degree arc.
 *
 ********************************************************************/
function globe_line(r, lon1, lat1, lon2, lat2)
{
    z1 = sin(lat1);
    x1 = cos(lat1) * cos(lon1);
    y1 = cos(lat1) * sin(lon1);
    
    z2 = sin(lat2);
    x2 = cos(lat2) * cos(lon2);
    y2 = cos(lat2) * sin(lon2);
    
    return globe_lxyz(r, x1, y1, z1, x2, y2, z2);
}