| 1026 | */ |
| 1027 | |
| 1028 | void makeRotationMatrix( double** B, const int _DIM, CBiteRnd& rrnd ) |
| 1029 | { |
| 1030 | double prod; |
| 1031 | int i, j, k; |
| 1032 | |
| 1033 | for( i = 0; i < _DIM; i++ ) |
| 1034 | { |
| 1035 | for( j = 0; j < _DIM; j++ ) |
| 1036 | { |
| 1037 | double unif = rrnd.get(); |
| 1038 | |
| 1039 | if( unif == 0.0 ) |
| 1040 | { |
| 1041 | unif = 1e-99; |
| 1042 | } |
| 1043 | |
| 1044 | double unif2 = rrnd.get(); |
| 1045 | |
| 1046 | if( unif2 == 0.0 ) |
| 1047 | { |
| 1048 | unif2 = 1e-99; |
| 1049 | } |
| 1050 | |
| 1051 | B[i][j] = sqrt(-2.0*log(unif)) * cos(2.0*M_PI*unif2); |
| 1052 | |
| 1053 | if( B[i][j] == 0 ) |
| 1054 | { |
| 1055 | B[i][j] = 1e-99; |
| 1056 | } |
| 1057 | } |
| 1058 | } |
| 1059 | |
| 1060 | for( i = 0; i < _DIM; i++ ) |
| 1061 | { |
| 1062 | for( j = 0; j < i; j++ ) |
| 1063 | { |
| 1064 | prod = 0.0; |
| 1065 | |
| 1066 | for( k = 0; k < _DIM; k++ ) |
| 1067 | { |
| 1068 | prod += B[k][i] * B[k][j]; |
| 1069 | } |
| 1070 | |
| 1071 | for( k = 0; k < _DIM; k++ ) |
| 1072 | { |
| 1073 | B[k][i] -= prod * B[k][j]; |
| 1074 | } |
| 1075 | } |
| 1076 | |
| 1077 | prod = 0.0; |
| 1078 | |
| 1079 | for( k = 0; k < _DIM; k++ ) |
| 1080 | { |
| 1081 | prod += B[k][i] * B[k][i]; |
| 1082 | } |
| 1083 | |
| 1084 | prod = sqrt( prod ); |
| 1085 | |