computes rotation matrices for real spherical harmonics
function (Rxyz, L) {
    
    var Nsh = (L+1)*(L+1);
    // allocate total rotation matrix
    var R = numeric.rep([Nsh,Nsh],0);
    
    // initialize zeroth and first band rotation matrices for recursion
    // Rxyz = [Rxx Rxy Rxz
    //         Ryx Ryy Ryz
    //         Rzx Rzy Rzz]
    //
    // zeroth-band (l=0) is invariant to rotation
    R[0][0] = 1;
    
    // the first band (l=1) is directly related to the rotation matrix
    var R_1 = numeric.rep([3,3],0);
    R_1[0][0] = Rxyz[1][1];
    R_1[0][1] = Rxyz[1][2];
    R_1[0][2] = Rxyz[1][0];
    R_1[1][0] = Rxyz[2][1];
    R_1[1][1] = Rxyz[2][2];
    R_1[1][2] = Rxyz[2][0];
    R_1[2][0] = Rxyz[0][1];
    R_1[2][1] = Rxyz[0][2];
    R_1[2][2] = Rxyz[0][0];
    
    R = numeric.setBlock(R, [1,1], [3,3], R_1);
    var R_lm1 = R_1;
    
    // compute rotation matrix of each subsequent band recursively
    var band_idx = 3;
    for (var l=2; l<L+1; l++) {
        
        var R_l = numeric.rep([(2*l+1),(2*l+1)],0);
        for (var m=-l; m<l+1; m++) {
            for (var n=-l; n<l+1; n++) {
                // compute u,v,w terms of Eq.8.1 (Table I)
                var d, denom, u, v, w;
                if (m==0) d = 1;
                else d = 0; // the delta function d_m0
                if (Math.abs(n)==l) denom = (2*l)*(2*l-1);
                else denom = (l*l-n*n);
                
                u = Math.sqrt((l*l-m*m)/denom);
                v = Math.sqrt((1+d)*(l+Math.abs(m)-1)*(l+Math.abs(m))/denom)*(1-2*d)*0.5;
                w = Math.sqrt((l-Math.abs(m)-1)*(l-Math.abs(m))/denom)*(1-d)*(-0.5);
                
                // computes Eq.8.1
                if (u!=0) u = u*U(l,m,n,R_1,R_lm1);
                if (v!=0) v = v*V(l,m,n,R_1,R_lm1);
                if (w!=0) w = w*W(l,m,n,R_1,R_lm1);
                R_l[m+l][n+l] = u + v + w;
            }
        }
        R = numeric.setBlock(R, [band_idx+1,band_idx+1], [band_idx+2*l+1,band_idx+2*l+1], R_l);
        R_lm1 = R_l;
        band_idx = band_idx + 2*l+1;
    }
    return R;
}