// Software    : Scilab >= 5.5
// Function    : mesharea()
// License     : CeCILL-C : http://www.cecill.info/licences
// Copyright(C): Samuel GOUGEON - France
// Version     : 1.0
// Date        : 2015-08-16
// Release     : http://fileexchange.scilab.org/toolboxes/370000
// History     :
//   2015-08-16 : First release on FileExchange
// 
//
// DISCLAIMER: This macro is provided as is, with NO WARRANTY of any kind.
//  Use it to your own risks, awareness, and responsability.

function ar = mesharea(X,Y,Z)
//
//Computes the area of any surface mapped with a {X,Y,Z} lattice of 4-vertices cells
//
// SYNTAX
// ------
//        mesharea  // Displays this help
// area = mesharea(X, Y, Z)
//
// PARAMETERS
// ----------
// X,Y,Z : 3 matrices of real decimal numbers, with same sizes: 
//         Cartesian coordinates of the nodes of the lattice mapping the surface.
//         {X(i,j), Y(i,j), Z(i,j)} locates the node (i,j).
//         Each node (not on a border or a corner) is linked to 4 neighboors.
//         Each cell of the lattice has 4 vertices set to nodes.
//         For a closed surface, the first and last rows of X|Y|Z must be equal,
//         as well as their first and last columns.
//         Nan values are accepted. Each node having at least one of its x|y|z
//         coordinates set to %nan is canceled. This allows for instance to make 
//         holes in the surface.
//         Folded surfaces are possible, when X, Y, and/or Z are not regular 
//         nor monotonic grids.
//
// area : single positive decimal number. The area of the mesh.
//        Each 4-vertices cell is split into 2 plan triangles. The area of the
//        cell is the sum of both triangles areas. The area of the mesh is the 
//        sum of all cells areas.
//        When a vertice is Nan, the area of the 4 triangles (half cells) 
//        touching it is ignored and not taken into account.
//
// REFERENCE
// ---------
//  Comments, scoring and bug reports are welcome on
//  http://fileexchange.scilab.org/toolboxes/370000
//
// EXAMPLE
// -------
//// Area of a sphere: Study of the error versus the surface sampling
//n = 4;
//for i = 1:8
//    long = linspace(0, 2*%pi, 2*n+1);     // sampling longitudes
//    lat  = linspace(-%pi/2, %pi/2, n+1)'; // sampling latitudes
//    X = cos(lat)*cos(long);
//    Y = cos(lat)*sin(long);
//    Z = sin(lat)*ones(long);
//    err(i) = mesharea(X,Y,Z)/(4*%pi) - 1; 
//    dAngle(i) = 180/n;  // angular step [°]
//    n = n*2;
//end
//clf, plot2d("ll", dAngle, -err)
//xtitle('mesharea() accuracy for a sphere', 'Angular sampling [°]', ..
//       'Relative error on sphere''s area')
//// => The relative error goes as the sampling step^2
//

    fname = "mesharea"
    if argn(2)==0
        head_comments(fname)
        ar = []
        return
    end
    
    // CHECKING PARAMETERS
    // -------------------
    if argn(2)~=3
        msg = _("%s: Wrong number of input arguments: %d expected.\n")
        error(msprintf(msg, fname, 3))
    end
    if type(X)~=1 | ~isreal(X) then
        msg = _("%s: Wrong type for input argument #%d: matrix of reals expected.\n")
        error(msprintf(msg, fname,1))
    end
    if type(Y)~=1 | ~isreal(Y) then
        msg = _("%s: Wrong type for input argument #%d: matrix of reals expected.\n")
        error(msprintf(msg, fname,2))
    end
    if type(Z)~=1 | ~isreal(Z) then
        msg = _("%s: Wrong type for input argument #%d: matrix of reals expected.\n")
        error(msprintf(msg, fname,3))
    end
    if or(size(Y)~=size(X)) then
        msg = _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n")
        error(msprintf(msg, fname,1,2))
    end
    if or(size(Z)~=size(X)) then
        msg = _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n")
        error(msprintf(msg, fname,1,3))
    end
    if size(X,1)<2 | size(X,2)<2
        msg = _("%s: Wrong size for input argument #1: Must be at least a (2x2) matrix.\n")
        error(msprintf(msg, fname,1))
    end

    // PROCESSING
    // ----------
    if X==[]
        ar = 0
        return
    end

    // First triangles
    dXv = X(2:$, 1:$-1) - X(1:$-1, 1:$-1)
    dYv = Y(2:$, 1:$-1) - Y(1:$-1, 1:$-1)
    dZv = Z(2:$, 1:$-1) - Z(1:$-1, 1:$-1)

    dXh = X(1:$-1, 2:$) - X(1:$-1, 1:$-1)
    dYh = Y(1:$-1, 2:$) - Y(1:$-1, 1:$-1)
    dZh = Z(1:$-1, 2:$) - Z(1:$-1, 1:$-1)
    
    dc = cross([dXh(:) dYh(:) dZh(:)]', [dXv(:) dYv(:) dZv(:)]')
    dn = sqrt(sum(dc.^2,"r"))
    i = find(~isnan(dn))
    ar = sum(dn(i))

    // 2nd triangles
    dXv = X(1:$-1, 2:$) - X(2:$, 2:$)
    dYv = Y(1:$-1, 2:$) - Y(2:$, 2:$)
    dZv = Z(1:$-1, 2:$) - Z(2:$, 2:$)

    dXh = X(2:$, 1:$-1) - X(2:$, 2:$)
    dYh = Y(2:$, 1:$-1) - Y(2:$, 2:$)
    dZh = Z(2:$, 1:$-1) - Z(2:$, 2:$)

    dc = cross([dXh(:) dYh(:) dZh(:)]', [dXv(:) dYv(:) dZv(:)]')
    dn = sqrt(sum(dc.^2,"r"))
    i = find(~isnan(dn))
    ar = ar + sum(dn(i))
    
    ar = ar/2
endfunction
