whole_func_string: 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;
  }
func_code_string: 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;
  }
func_documentation_string: computes rotation matrices for real spherical harmonics