```
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <ctime>
#include <cmath>
using namespace std;
const float M_2PIf = 6.283185307179586476925286766559f;
const float M_PIf = 3.141592653589793238462643383279f;
struct vector {
float x, y, z;
};
inline float function_vec(const vector &left, const vector &right)
{
return(left.x * right.x + left.y * right.y + left.z * right.z);
}
class Matrix3x3 {
private:
float tomb[3][3]; // 3x3 matrix
public:
void init(vector &a){
tomb[0][0] = a.x * a.x;
tomb[0][1] = a.x * a.y;
tomb[0][2] = a.x * a.z;
tomb[1][0] = a.y * a.x;
tomb[1][1] = a.y * a.y;
tomb[1][2] = a.y * a.z;
tomb[2][0] = a.z * a.x;
tomb[2][1] = a.z * a.y;
tomb[2][2] = a.z * a.z;
}
Matrix3x3 &operator +=(Matrix3x3 b){
int k,l;
for(k=0; k<3; k++)
for(l=0; l<3; l++)
tomb[k][l] += b.tomb[k][l];
return *this;
}
Matrix3x3 &operator/=(float c){
int k,l;
for(k=0; k<3; k++)
for(l=0; l<3; l++)
tomb[k][l] /= c;
return *this;
}
float element(int i,int j) const {
return tomb[i][j];
}
};
ostream& operator << (ostream& output, const Matrix3x3& s){
output << "[\n";
for(int i=0; i<3; i++){
for(int j=0; j<3; j++)
output << "\t" << s.element(i,j) << "\t";
output << "\n";
}
output << "]";
return output;
}
// Function prototypes
float rand(float min, float max);
float vec_angle(const vector &left, const vector &right);
vector function(float m[][3], vector v);
// main function
int main(void) {
int N = 100;
float epsilon = 1.0;
float cosangle[N];
// >>>>>>>>> HISTOGRAMM <<<<<<<<<<<<
double hmin = 0.0;
double hmax = 1.0;
int hbins = 10;
double* h;
h = new double[hbins];
for ( int bin = 0 ; bin<hbins ; bin++ ) h[bin] = 0.0;
// >>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<
Matrix3x3 m[N];
vector vec = {0.0,0.0,1.0};
vector vec0 = {0.0,0.0,1.0};
int i, j;
srand((unsigned)time(NULL));
ofstream ofile;
ofile.open("datavecA.dat");
for (i = 0; i < N; i++)
{
// >>>>>>>>> ________________ <<<<<<<<<<<
float mX[3][3] = {
{1.0f, 0.0f, 0.0f},
{0.0f, 1.0f, 0.0f},
{0.0f, 0.0f, 1.0f}
};
float mY[3][3] = {
{1.0f, 0.0f, 0.0f},
{0.0f, 1.0f, 0.0f},
{0.0f, 0.0f, 1.0f}
};
float mZ[3][3] = {
{1.0f, 0.0f, 0.0f},
{0.0f, 1.0f, 0.0f},
{0.0f, 0.0f, 1.0f}
};
// >>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<
switch (rand()%3)
{
case 0:
vec = function(mX, vec);
break;
case 1:
vec = function(mY, vec);
break;
case 2:
vec = function(mZ, vec);
break;
}
cosangle[i] = vec_angle(vec0,vec);
ofile << cosangle[i] << endl;
int bin=(int) ((cosangle[i]-hmin)/(hmax-hmin)*hbins);
if ( bin<0 ) bin=0;
if ( bin>hbins-1 ) bin=hbins-1;
h[bin]++;
m[i].init(vec);
}
ofile.close();
ofstream file;
file.open("datavecB.dat");
file << "#position\tvalue\t\thalf-binwidth\terror" << endl;
for ( int bin=0 ; bin<hbins ; bin++ )
{
double X=hmin+(bin+0.5)*(hmax-hmin)/hbins;
double dX=0.5*(hmax-hmin)/hbins;
double Y=h[bin];
double dY=sqrt(h[bin]);
file << " " << X << "\t\t" << Y << "\t\t" << dX << "\t\t" << dY << endl;
}
file.close();
delete[] h;
/****************************************************/
Matrix3x3 sum = m[0];
for (i=0;i<N;i++) {
sum += m[i];
}
sum /= N;
cout << sum << endl;
return 0;
}
float rand(float min, float max)
{
return min + (max - min) * rand() / (float)RAND_MAX;
}
vector function(float m[][3], vector v)
{
vector returnv;
returnv.x = m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z;
returnv.y = m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z;
returnv.z = m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z;
return returnv;
}
float vec_angle(const vector &left, const vector &right)
{
vector v1 = left, v2 = right;
float dot = function_vec(v1,v2);
if(fabsf(dot)>=1.0f) dot = 0.999f;
return acosf(dot);
}
```