#################################################################################################
#
#  3 次元の orthonormal matrix を、fine grid で作る。
#  
#  注意:
#  ※ これらで作られる matrix の示すベクトルは一様には分布していないことに注意。
#  ※ 軸の選択についての不定性(鏡像・軸の交換)について考えなければ、基本的にこれでよさそう。
#  ※ 3 次元で、30(31) 分割なら、29791 個の行列となり、R のオブジェクトのサイズも 3Mb 程度でよい。
#
#  パラメタ:
#  ※ 回転角の中心値( th0, ph0, lm0 ) を指定できる。
#  ※ 計算する範囲( range ) を radian で指定できる。
#  ※ メッシュの間隔を個数( nmesh )で指定できる。
#
#################################################################################################

MakeFineGridOrthonormalMatrix <- function(dim=3, nmesh=30, th0=pi/4, ph0=pi/4, lm0=pi/4, range=pi/2){

  matrixlist <- list(NULL)
  angles     <- list(NULL)
  counter    <- 0

# 2D
  if(       identical( dim, 2 ) ){
    for( theta in seq( th0 - range/2, th0 + range/2, range/nmesh ) ){
      W <- matrix( c( cos(theta), -sin(theta), sin(theta),cos(theta)),2 )
      counter <- counter + 1
      matrixlist[[counter]] <- W
      angles[[counter]]     <- theta
    }

# 3D
  }else if( identical( dim, 3 ) ){
    for( theta in seq(    th0 - range/2, th0 + range/2, range/nmesh ) ){
      for( phi in seq(    ph0 - range/2, ph0 + range/2, range/nmesh ) ){ 
        for( lamd in seq( lm0 - range/2, lm0 + range/2, range/nmesh ) ){
          W <- 
            matrix( c( cos(lamd), sin(lamd), 0, -sin(lamd), cos(lamd), 0, 0, 0, 1 ), 3)       %*% # z-axis
              matrix( c( 1, 0, 0, 0, cos(theta), sin(theta), 0, -sin(theta), cos(theta) ), 3) %*% # x-axis
                matrix( c( sin(phi), 0, cos(phi), 0, 1, 0,  cos(phi), 0, -sin(phi) ), 3)          # y-axis
          counter <- counter + 1
          matrixlist[[counter]] <- W
          angles[[counter]]     <- c( theta, phi, lamd )
        }
      }
    }

# その他
  }else{
    print("This function supports 2D and 3D only\n")
    return(0)
  }

  return(list(matrix=matrixlist,angles=angles))
}