/////////////////////////////////////////
//Matrix.h
#ifndef _MATRIX_H_NOU_
#define _MATRIX_H_NOU_

#include <string.h>
#include <math.h>
#include <windows.h>
#include <gl/gl.h>

#ifndef RAD
#define RAD (180/3.14159265358979323846)
#endif

class Vector
{
private:
    float vector[3];
public:
    Vector(){vector[0]=vector[1]=vector[2]=0;}
    float operator[](const int i)const{return vector[i];}
    void operator()(int i,float x){vector[i]=x;}
    void operator()(float *i){vector[0]=i[0];vector[1]=i[1];vector[2]=i[2];}
    void Set(float *i){vector[0]=i[0];vector[1]=i[1];vector[2]=i[2];}
    Vector operator+(const Vector &vec);//sucet dvoch vektorov
    Vector operator-(const Vector &vec);//rozdiel vektorov
    Vector operator*(const Vector &vec);//vektorovy sucin vektorov
    float Skalar(const Vector &vec);//skalarny sucin
    void Normalize();//znormalizuje vektor
    float UholR(Vector &vec);//vrati uhol v Radianoch
    float UholS(Vector &vec)//vrati uhol dvoch vektorov v stupnoch
    {
        return this->UholR(vec)*RAD;
    }
    float Velkost(){return sqrt(vector[0]*vector[0]+vector[1]*vector[1]+vector[2]*vector[2]);}
};

class Bod
{
private:
    float bod[3];
public:
    Bod(){bod[0]=bod[1]=bod[2]=0;}
    float operator[](int i)const{return bod[i];}
    void operator()(float *i){bod[0]=i[0];bod[1]=i[1];bod[2]=i[2];}
    void Set(float *i){bod[0]=i[0];bod[1]=i[1];bod[2]=i[2];}
    void operator()(int i, float x){if(i>=0 && i<=2)bod[i]=x;}
    Vector operator-(const Bod &bod2)//rozdiel dvoch bodov
    {
        Vector v;
        v(0 ,bod[0]-bod2[0]);
        v(1 ,bod[1]-bod2[1]);
        v(2 ,bod[2]-bod2[2]);
        return v;
    }
    float Delta(Bod &bod)//vzdialenost dvoch bodov
    {
        Vector rozdiel;
        rozdiel=*this-bod;
        return (float)sqrt(rozdiel[0]*rozdiel[0]+rozdiel[1]*rozdiel[1]+rozdiel[2]*rozdiel[2]);
    }
};

class Matrix
{
private:
    float matica[16];
public:
    Matrix();
    Matrix operator*(const Matrix &m2);//nasobenie matic
    Vector operator*(const Vector &vec);//nasobenie vectora a matice
    Bod operator*(const Bod &bod);//nasobenie bodu a matice
    Matrix& operator*=(Matrix &m2){*this=m2**this;return *this;}
    float operator[](int i)const{return matica[i];}
    void operator()(int i,float x){if(i>=0 && i<=15)matica[i]=x;}
    void Inverse();//zinvertizuje maticu
    float Determinant();//pocita determinant matice
    void Identity()//zmeni na jednotkovu maticu
    {
        memset(matica,0,sizeof(float)*16);
        matica[0]=matica[5]=matica[10]=matica[15]=1;
    }
    void operator()(float *m){memcpy(matica,m,sizeof(float)*16);}//nastavime z vektora hodnotu matice
    void Rotate(float x,float y,float z);//Eulerova rotacia
    void Translate(float x,float y,float z)
    {
        matica[12]=x;  matica[13]=y;    matica[14]=z;
    }
    void LoadGL(){glLoadMatrixf(matica);}//nacita maticu do OGL
    void GetGL(GLenum type=GL_MODELVIEW_MATRIX){glGetFloatv(type,matica);}
    void GetMatrix(float *mat){memcpy(mat,matica,sizeof(float)*16);}
};


#endif

/////////////////////////////////////////
//Matrix.cpp
#include "Matrix.h"
#include <stdio.h>

Vector Vector::operator+(const Vector &vec)//sucet dvoch vektorov
{
    Vector v;
    v(0, vec[0]+vector[0]);
    v(1, vec[1]+vector[1]);
    v(2, vec[2]+vector[2]);
    return v;
}

Vector Vector::operator-(const Vector &vec)//rozdiel vektorov
{
    Vector v;
    v(0, vector[0]-vec[0]);
    v(1, vector[1]-vec[1]);
    v(2, vector[2]-vec[2]);
    return v;
}

Vector Vector::operator*(const Vector &vec)//vektorovy sucin vektorov
{
    Vector v;
    v(0, vector[1]*vec[2]-vector[2]*vec[1]);
    v(1, vector[2]*vec[0]-vector[0]*vec[2]);
    v(2, vector[0]*vec[1]-vector[1]*vec[0]);
    return v;
}

float Vector::Skalar(const Vector &vec)
{
    return vector[0]*vec[0]+vector[1]*vec[1]+vector[2]*vec[2];
}

void Vector::Normalize()//znormalizuje vektor
{
    float d=(float)sqrt(vector[0]*vector[0]+vector[1]*vector[1]+vector[2]*vector[2]);
    vector[0]/=d;
    vector[1]/=d;
    vector[2]/=d;
}

float Vector::UholR(Vector &vec)//vrati uhol v Radianoch
{
    Vector u,v;
    u=*this;
    v=vec;
    u.Normalize();
    v.Normalize();
    return (float)acos(u.Skalar(v));
}

Matrix::Matrix()
{
    memset(matica,0,sizeof(float)*16);
    matica[0]=matica[5]=matica[10]=matica[15]=1;
}

Matrix Matrix::operator*(const Matrix &m2)//nasobenie matic
{
    Matrix mt;
    float *m1=matica;
    mt(0, m1[0]*m2[0]+m1[1]*m2[4]+m1[2]*m2[8]+m1[3]*m2[12]);
    mt(1, m1[0]*m2[1]+m1[1]*m2[5]+m1[2]*m2[9]+m1[3]*m2[13]);
    mt(2, m1[0]*m2[2]+m1[1]*m2[6]+m1[2]*m2[10]+m1[3]*m2[14]);
    mt(3, m1[0]*m2[3]+m1[1]*m2[7]+m1[2]*m2[11]+m1[3]*m2[15]);

    mt(4, m1[4]*m2[0]+m1[5]*m2[4]+m1[6]*m2[8]+m1[7]*m2[12]);
    mt(5, m1[4]*m2[1]+m1[5]*m2[5]+m1[6]*m2[9]+m1[7]*m2[13]);
    mt(6, m1[4]*m2[2]+m1[5]*m2[6]+m1[6]*m2[10]+m1[7]*m2[14]);
    mt(7, m1[4]*m2[3]+m1[5]*m2[7]+m1[6]*m2[11]+m1[7]*m2[15]);

    mt( 8, m1[8]*m2[0]+m1[9]*m2[4]+m1[10]*m2[8]+m1[11]*m2[12]);
    mt( 9, m1[8]*m2[1]+m1[9]*m2[5]+m1[10]*m2[9]+m1[11]*m2[13]);
    mt(10, m1[8]*m2[2]+m1[9]*m2[6]+m1[10]*m2[10]+m1[11]*m2[14]);
    mt(11, m1[8]*m2[3]+m1[9]*m2[7]+m1[10]*m2[11]+m1[11]*m2[15]);

    mt(12, m1[12]*m2[0]+m1[13]*m2[4]+m1[14]*m2[8]+m1[15]*m2[12]);
    mt(13, m1[12]*m2[1]+m1[13]*m2[5]+m1[14]*m2[9]+m1[15]*m2[13]);
    mt(14, m1[12]*m2[2]+m1[13]*m2[6]+m1[14]*m2[10]+m1[15]*m2[14]);
    mt(15, m1[12]*m2[3]+m1[13]*m2[7]+m1[14]*m2[11]+m1[15]*m2[15]);

    return mt;
}

void Matrix::Inverse()//inverzna matica
{
    float minor[16];
    float *m=matica;
    //vyratame minor zaroven cofaktor aj otocenie po uhlopriecke
    minor[ 0]=m[5]*m[10]*m[15]+m[6]*m[11]*m[13]+m[7]*m[9]*m[14]-m[7]*m[10]*m[13]-m[6]*m[9]*m[15]-m[5]*m[11]*m[14];
    minor[ 4]=-(m[4]*m[10]*m[15]+m[6]*m[11]*m[12]+m[7]*m[8]*m[14]-m[7]*m[10]*m[12]-m[6]*m[8]*m[15]-m[4]*m[11]*m[14]);
    minor[ 8]=m[4]*m[9]*m[15]+m[5]*m[11]*m[12]+m[7]*m[8]*m[13]-m[7]*m[9]*m[12]-m[5]*m[8]*m[15]-m[4]*m[11]*m[13];
    minor[12]=-(m[4]*m[9]*m[14]+m[5]*m[10]*m[12]+m[6]*m[8]*m[13]-m[6]*m[9]*m[12]-m[5]*m[8]*m[14]-m[4]*m[10]*m[13]);

    minor[ 1]=-(m[1]*m[10]*m[15]+m[2]*m[11]*m[13]+m[3]*m[9]*m[14]-m[3]*m[10]*m[13]-m[2]*m[9]*m[15]-m[1]*m[11]*m[14]);
    minor[ 5]=m[0]*m[10]*m[15]+m[2]*m[11]*m[12]+m[3]*m[8]*m[14]-m[3]*m[10]*m[12]-m[2]*m[8]*m[15]-m[0]*m[11]*m[14];
    minor[ 9]=-(m[0]*m[9]*m[15]+m[1]*m[11]*m[12]+m[3]*m[8]*m[13]-m[3]*m[9]*m[12]-m[1]*m[8]*m[15]-m[0]*m[11]*m[13]);
    minor[13]=m[0]*m[9]*m[14]+m[1]*m[10]*m[12]+m[2]*m[8]*m[13]-m[2]*m[9]*m[12]-m[1]*m[8]*m[14]-m[0]*m[10]*m[13];

    minor[ 2]=m[1]*m[6]*m[15]+m[2]*m[7]*m[13]+m[3]*m[5]*m[14]-m[3]*m[6]*m[13]-m[2]*m[5]*m[15]-m[1]*m[7]*m[14];
    minor[ 6]=-(m[0]*m[6]*m[15]+m[2]*m[7]*m[12]+m[3]*m[4]*m[14]-m[3]*m[6]*m[12]-m[2]*m[4]*m[15]-m[0]*m[7]*m[14]);
    minor[10]=m[0]*m[5]*m[15]+m[1]*m[7]*m[12]+m[3]*m[4]*m[13]-m[3]*m[5]*m[12]-m[1]*m[4]*m[15]-m[0]*m[7]*m[13];
    minor[14]=-(m[0]*m[5]*m[14]+m[1]*m[6]*m[12]+m[2]*m[4]*m[13]-m[2]*m[5]*m[12]-m[1]*m[4]*m[14]-m[0]*m[6]*m[13]);

    minor[ 3]=-(m[1]*m[6]*m[11]+m[2]*m[7]*m[9]+m[3]*m[5]*m[10]-m[3]*m[6]*m[9]-m[2]*m[5]*m[11]-m[1]*m[7]*m[10]);
    minor[ 7]=m[0]*m[6]*m[11]+m[2]*m[7]*m[8]+m[3]*m[4]*m[10]-m[3]*m[6]*m[8]-m[2]*m[4]*m[11]-m[0]*m[7]*m[10];
    minor[11]=-(m[0]*m[5]*m[11]+m[1]*m[7]*m[8]+m[3]*m[4]*m[9]-m[3]*m[5]*m[8]-m[1]*m[4]*m[11]-m[0]*m[7]*m[9]);
    minor[15]=m[0]*m[5]*m[10]+m[1]*m[6]*m[8]+m[2]*m[4]*m[9]-m[2]*m[5]*m[8]-m[1]*m[4]*m[10]-m[0]*m[6]*m[9];

    float det=Determinant();
    //vydelit adjoint / |A|
    if(det)//proti deleniu nulou
    {
        for(int i=0;i<16;i++)
            matica[i]=minor[i]/det;
    }
}

float Matrix::Determinant()
{
    //45 operacii najprv vyratame 2x2 determinanty
    float a1=matica[10]*matica[15]-matica[11]*matica[14];
    float a2=matica[ 9]*matica[15]-matica[11]*matica[13];
    float a3=matica[ 9]*matica[14]-matica[10]*matica[13];
    float a4=matica[ 7]*matica[ 2]-matica[ 6]*matica[ 3];
    float a5=matica[ 7]*matica[ 1]-matica[ 5]*matica[ 3];
    float a6=matica[ 6]*matica[ 1]-matica[ 5]*matica[ 2];

    //vratime cely determinant
    return matica[ 0]*(matica[ 5]*a1-matica[ 6]*a2+matica[ 7]*a3)-
           matica[ 4]*(matica[ 1]*a1-matica[ 2]*a2+matica[ 3]*a3)+
           matica[ 8]*(matica[13]*a4-matica[14]*a5+matica[15]*a6)-
           matica[12]*(matica[ 9]*a4-matica[10]*a5+matica[11]*a6);
}

Vector Matrix::operator*(const Vector &vec)//nasobenie vectora a matice
{
    Vector v;
    v(0, vec[0] * matica[0] + vec[1] * matica[4] + vec[2] * matica[8]);
    v(1, vec[0] * matica[1] + vec[1] * matica[5] + vec[2] * matica[9]);
    v(2, vec[0] * matica[2] + vec[1] * matica[6] + vec[2] * matica[10]);
    return v;
}

Bod Matrix::operator*(const Bod &bod)//nasobenie bodu a matice
{
    Bod b;
    b(0, bod[0] * matica[0] + bod[1] * matica[4] + bod[2] * matica[8] + matica[12]);
    b(1, bod[0] * matica[1] + bod[1] * matica[5] + bod[2] * matica[9] + matica[13]);
    b(2, bod[0] * matica[2] + bod[1] * matica[6] + bod[2] * matica[10] + matica[14]);
    return b;
}

void Matrix::Rotate(float x, float y, float z)// Hlavní funkce pro rotaci
{
    glPushMatrix();
    glLoadMatrixf(matica);
    if(z)
        glRotatef(z*RAD,0,0,1);
    if(y)
        glRotatef(y*RAD,0,1,0);
    if(x)
        glRotatef(x*RAD,1,0,0);
    glGetFloatv(GL_MODELVIEW_MATRIX,matica);
    glPopMatrix();
}