/********************************************************************
 * globe plot class
 *
 * Example:
    import globe_class;
    G = new Globe;
    G->surface();
    G->gshhs("data/gshhs_c.b");
    G->grid(20);
    G->show();
 *
 ********************************************************************/

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

D2R = 0.017453292519943;    // degree to radius

SCALE = 0.35;               // globe to image ratio

class Globe {

    c_width = 600;               // image size
    c_height = 600;
    
    c_rotatex = -60;            // Atlentic centered
    c_rotatez = -90;
        
    c_radius = SCALE * c_width;   // globe size

    c_render = zegraph("render");
    c_scene = zegraph("scene");
    c_root = zegraph("node");
    
    c_render:size(c_width, c_height);
    c_scene:viewport(0, 0, c_width, c_height);
    c_render:add(c_scene);
    c_scene:root(c_root);
    
    c_light = zegraph("light");
    c_light:position(0, 0.25, 1);
    c_light:ambient(0.5, 0.5, 0.5);
    
    c_node = zegraph("node");
    c_node:rotatez(c_rotatez);
    c_node:rotatex(c_rotatex);
    
    c_root:add(c_light, c_node);

    c_font = zegraph("font");
    
   /***************************************************************
    * set size
    ***************************************************************/
    function size(width, height)
    {
        global:c_width = width;
        global:c_height = height;
        c_render:size(width, height);
        c_scene:viewport(0, 0, width, height);
        global:c_radius = SCALE * width;
    }  
    
   /***************************************************************
    * add object to the node
    ***************************************************************/
    function add(shape, transform=true)
    {
        if (transform)
            c_node:add(shape);
        else
            c_root:add(shape);
    }  
    
   /***************************************************************
    * rotation
    ***************************************************************/
    function rotate(rotatex, rotatez)
    {
        global:c_rotatex = rotatex;
        global:c_rotatez = rotatez;
        c_node:rotatez(rotatez);
        c_node:rotatex(rotatex);
    }  
    
   /****************************************************************
    * move globbe up-down or left-right
    ****************************************************************/   
    function move(dx, dy)
    {
        c_root:translate(dx, dy, 0);
    }
    
   /***************************************************************
    * rotation
    ***************************************************************/
    function view(lat, lon)
    {
        global:c_rotatex = -(90 - lat);
        global:c_rotatez = -(90 + lon);
        c_node:rotatez(c_rotatez);
        c_node:rotatex(c_rotatex);
    }  
    
   /***************************************************************
    * show the globe
    ***************************************************************/
    function show()
    {
        window(c_width, c_height, "show_callback");

        function show_callback(hwnd, msg, wp, hwp, lwp, lp, hlp, llp)
        {
            if (msg == "PAINT") {
                c_render:towindow(hwnd);
            }
            else if (msg == "SIZE") {
                hide();
            }
        }
    }

   /***************************************************************
    * save image to file
    ***************************************************************/
    function save(fname)
    {
        c_render:tofile(fname);
    }
    
   /****************************************************************
    * globe surface
    ****************************************************************/
    function surface(iteration=5)
    {
        polygon = zegraph("polygon");
        polygon:color(0, 0, .5);
        c_node:add(polygon);
        polygon:sphere(c_radius, iteration);
        return polygon;
    }
    
   /****************************************************************
    * create a globe using image for texutre 
    ****************************************************************/
    function texture(fname)
    {
        node = zegraph("node");
        texture = zegraph("texture");
        polygon = zegraph("polygon");
        polygon:color(1, 1, 1);
        node:add(texture, polygon);
        node:closed(true);
        c_node:add(node);
        texture:image(fname, -1, -1, -1);
        polygon:sphere(c_radius, 5, true);
    }
    
   /***************************************************************
    * draw grid lines
    ****************************************************************/
    function grid(deg, width=1)
    {
        deg = integer(deg);
        assert(deg > 1);
        r = 1.001*c_radius;
        
        node = zegraph("node");
        node:color(.5, .5, .5);
        c_node:add(node);
        
        n = (180 - 2*deg)/2 + 1;
        
        D = matrix("double", n, 1);
        D:fill(-90+deg, 2);
        
        XYZ = matrix("double", n, 3);
        XYZ[null,1] = D;
        XYZ[null,2] = 0;
    
        for (lon = 0; lon < 360; lon += deg) {
            XYZ[null,0] = lon;
            line = zegraph("line");
            node:add(line);
            line:type("strip");
            line:solid(width);
            xyz = zegraph("vertex");
            line:vertex(xyz);
            ptr = XYZ:ptr();
            xyz:add(ptr[0], ptr[1]);
            xyz:xy2sphere(r);
            nor = xyz:normal(1);
            line:normal(nor);
        }
        
        n = 181;
        
        D:resize(n, 1);
        D:fill(0, 2);
        
        XYZ:resize(n, 3);
        XYZ[null,0] = D;
        XYZ[null,2] = 0;
            
        for (lat = -90+deg; lat <= 90-deg; lat += deg) {
            XYZ[null,1] = lat;
            line = zegraph("line");
            node:add(line);
            line:type("strip");
            line:solid(width);
            xyz = zegraph("vertex");
            line:vertex(xyz);
            ptr = XYZ:ptr();
            xyz:add(ptr[0], ptr[1]);
            xyz:xy2sphere(r);
            nor = xyz:normal(1);
            line:normal(nor);
        }
    
        return node;
    }
    
   /****************************************************************
    * make coastlines using GSHHS data
    ****************************************************************/
    function gshhs(fname)
    {
        r = 1.001*c_radius;
        
        line = zegraph("line");
        c_node:add(line);
        xyz = zegraph("vertex");
        nor = zegraph("vertex");
        line:type("lines");
        line:color(.8, .8, .8);
        line:vertex(xyz);
        
        D = matrix("double");
        D:readgshhs(fname, "land");
        D:insert(2, 0);
        ptr = D:ptr();
        xyz:add(ptr[0], ptr[1]);
        xyz:xy2sphere(r);
        nor = xyz:normal(1);
        line:normal(nor);
        
        return line;
    }
    
   /****************************************************************
    * draw a line on surface.
    ****************************************************************/
    function line(lon1, lat1, lon2, lat2)
    {
        lon1 *= D2R;
        lon2 *= D2R;
        lat1 *= D2R;
        lat2 *= D2R;

        z1 = sin(lat1);
        r = cos(lat1);
        x1 = r * cos(lon1);
        y1 = r * sin(lon1);
        
        z2 = sin(lat2);
        r = cos(lat2);
        x2 = r * cos(lon2);
        y2 = r * sin(lon2);

        r = 1.001*c_radius;
        
        line = zegraph("line");
        c_node:add(line);
        vertex = zegraph("vertex");
        normal = zegraph("vertex");
        line:type("strip");
        line:vertex(vertex);
        line:normal(normal);
        
        dx = real(x2 - x1);
        dy = real(y2 - y1);
        dz = real(z2 - z1);
        d = sqrt(dx*dx + dy*dy + dz*dz);
        assert(d > 0);
        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;
    }
    
   /****************************************************************
    * text layout
    ****************************************************************/ 
    function text(x, y, string, size=12, halign=0, valign=-1)
    {
        text = zegraph("text");
        text:font(c_font, size);
        text:string(string);
        size = text:size();
        h = -size[0]/2;
        v = -size[1]/2;
        if (halign < 0) h = 0;
        if (halign > 0) h = -size[0];
        if (valign < 0) v = 0;
        if (valign > 0) v = -size[1];
        text:layout(x+h, y+v, 0,  x+h+1, y+v, 0,  x+h, y+v+1, 0);
        c_root:add(text);
        return text;
    }
    
   /****************************************************************
    * lon-lat to x-y-z projection
    ****************************************************************/
    function projection(lon, lat)
    {
        lon *= D2R;
        lat *= D2R;
        z = c_radius * sin(lat);
        r = c_radius * cos(lat);
        x = r * cos(lon);
        y = r * sin(lon);
        return [x, y, z];
    }
    
    /****************************************************************
     * greate circle distance on unit globe
     ****************************************************************/
    function distance(lon1, lat1, lon2, lat2)
    {
        lon1 *= D2R;
        lat1 *= D2R;
        lon2 *= D2G;
        lat2 *= D2R;
        return acos(cos(lat1)*cos(lat2)*cos(lon1-lon2) + sin(lat1)*sin(lat2));
    }
    
   /****************************************************************
    * animate the globe and control rotation by arrow kyes and mouse
    ********************************************************************/
    function animate()
    {
        rotx0 = c_rotatex;
        rotz0 = c_rotatez;
        rotx = rotx0;
        rotz = rotz0;
        
        scx = c_width / 2;
        scy = c_height / 2;
    
        x0 = 0;
        y0 = 0;
        cx = 0;
        cy = 0;
        drag = 0;
        deg = 15;
        clock = null;
        key = 39;
        
        window(c_width, c_height, "animate_callback");
    
        function animate_callback(hwnd, msg, wp, hwp, lwp, lp, hlp, llp)
        {
            if (msg == "PAINT") {
                c_render:towindow(hwnd);
            }
            else if (msg =="KEYUP") {
                if (wp >= 37 && wp <= 40) {
                    key = wp;
                    if (wp == 37)               // left arrow
                        global:rotz -= deg;
                    else if (wp == 38)          // up arrow
                        global:rotx -= deg;
                    else if (wp == 39)          // right arrow
                        global:rotz += deg;
                    else                        // down arrow
                        global:rotx += deg;
                    
                    c_node:rotatez(rotz);
                    c_node:rotatex(rotx);
                    c_render:towindow(hwnd);
                }
            }
            else if (msg == "LBUTTONDOWN") {
                global:x0 = hlp;
                global:y0 = llp;
                global:drag = 1;
            }
            else if (msg == "LBUTTONUP") {
                global:drag = 0;
            }
            else if (msg == "RBUTTONUP") {
                global:drag = 0;
                c_node:rotatez(rotz0);
                c_node:rotatex(rotx0);
                c_render:towindow(hwnd);
            }
            else if (msg == "MOUSEMOVE") {
                if (drag > 0) {
                    if (cx != x0 && cx != hlp) {
                        if ((x0 - hlp) / (cx - hlp) < 0) {
                            global:x0 = hlp;
                            global:cx = x0;
                            global:y0 = llp;
                            global:cy = y0;
                            return;
                        }
                    }
                    if (cy != y0 && cy != llp) {
                        if ((y0 - llp) / (cy - llp) < 0) {
                            global:x0 = hlp;
                            global:cx = x0;
                            global:y0 = llp;
                            global:cy = y0;
                            return;
                        }
                    }
                    a = arcball(x0, y0, hlp, llp, scx, scy);
                    global:rotx = rotx + a[0];
                    global:rotz = rotz + a[1];
                    c_node:rotatez(rotz);
                    c_node:rotatex(rotx);
                    c_render:towindow(hwnd);
                    global:cx = hlp;
                    global:cy = llp;
                }
            }
            else if (msg =="LBUTTONDBLCLK") {
                if (clock == null) {
                    timer(30);
                    global:clock = true;
                }
                else {
                    timer(0);
                    global:clock = null;
                }
            }
            else if (msg == "TIMER") {
                if (key == 37)
                    global:rotz -= 2;
                else if (key == 38)
                    global:rotx -= 2;
                else if (key == 39)
                    global:rotz += 2;
                else
                    global:rotx += 2;
                    
                c_node:rotatez(rotz);
                c_node:rotatex(rotx);
                c_render:towindow(hwnd);
            }
            else if (msg == "SIZE") {
                hide();
            }
        }
    }
    
   /********************************************************************
    * Arcball emulation
    ********************************************************************/
    function arcball(x0, y0, x1, y1, cx, cy)
    {
        x0 -= cx;
        x0 /= c_radius;
        y0 -= cy;
        y0 /= -c_radius;
        a = x0*x0 + y0*y0;
    
        if (a > 1.0) {
            a = sqrt(a);
            x0 /= a;
            y0 /= a;
            a = 1.0;
        }
        z0 = sqrt(1 - a);
    
        x1 -= cx;
        x1 /= c_radius;
        y1 -= cy;
        y1 /= -c_radius;
        a = x1*x1 + y1*y1;
        if (a > 1.0) {
            a = sqrt(a);
            x1 /= a;
            y1 /= a;
            a = 1.0;
        }
        z1 = sqrt(1 - a);
        
        a = direction(0, 0, 0, z0, x0, y0);
        b = direction(0, 0, 0, z1, x1, y1);
    
        return [a[0] - b[0], b[1] - a[1]];
    }
}