nl_tokens: computes rotation matrices for real spherical harmonics
pl_tokens: function ( Rxyz , L ) { var Nsh = ( L + 1 ) * ( L + 1 ) ; // allocate total rotation matrix var R = numeric . <mask> ( [ 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 ; }