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;
}
function (Rxyz, L) {
    
    var Nsh = (L+1)*(L+1);
    
    var R = numeric.rep([Nsh,Nsh],0);
    
    
    
    
    
    
    
    R[0][0] = 1;
    
    
    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;
    
    
    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++) {
                
                var d, denom, u, v, w;
                if (m==0) d = 1;
                else d = 0; 
                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);
                
                
                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;
}
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 ; }
computes rotation matrices for real spherical harmonics
computes rotation matrices for real spherical harmonics