computeRealSH computes real spherical harmonics up to order N
function (N, data) {
    
    var azi = new Array(data.length);
    var elev = new Array(data.length);
    
    for (var i=0; i<data.length; i++) {
        azi[i] = data[i][0];
        elev[i] = data[i][1];
    }
    
    var factorials = new Array(2*N+1);
    var Ndirs = azi.length;
    var Nsh = (N+1)*(N+1);
    var leg_n_minus1 = 0;
    var leg_n_minus2 = 0;
    var leg_n;
    var sinel = numeric.sin(elev);
    var index_n = 0;
    var Y_N = new Array(Nsh);
    var Nn0, Nnm;
    var cosmazi, sinmazi;
    
    // precompute factorials
    for (var i = 0; i < 2*N+1; i++) factorials[i] = factorial(i);
    
    for (var n = 0; n<N+1; n++) {
        if (n==0) {
            var temp0 = new Array(azi.length);
            temp0.fill(1);
            Y_N[n] = temp0;
            index_n = 1;
        }
        else {
            leg_n = recurseLegendrePoly(n, sinel, leg_n_minus1, leg_n_minus2);
            Nn0 = Math.sqrt(2*n+1);
            for (var m = 0; m<n+1; m++) {
                if (m==0) Y_N[index_n+n] = numeric.mul(Nn0,leg_n[m]);
                else {
                    Nnm = Nn0*Math.sqrt( 2 * factorials[n-m]/factorials[n+m] );
                    cosmazi = numeric.cos(numeric.mul(m,azi));
                    sinmazi = numeric.sin(numeric.mul(m,azi));
                    Y_N[index_n+n-m] = numeric.mul(Nnm, numeric.mul(leg_n[m], sinmazi));
                    Y_N[index_n+n+m] = numeric.mul(Nnm, numeric.mul(leg_n[m], cosmazi));
                }
            }
            index_n = index_n+2*n+1;
        }
        leg_n_minus2 = leg_n_minus1;
        leg_n_minus1 = leg_n;
    }
    
    return Y_N;
}