#ifndef __VECTOR_H
#define __VECTOR_H

struct VECTOR {
    float x, y, z;

    VECTOR() {}
    VECTOR(float _x, float _y, float _z) {
        x = _x; y = _y; z = _z;
    }

    inline operator +=(const VECTOR &v);
    inline operator -=(const VECTOR &v);
    inline operator *=(float f);
    inline operator /=(float f);
};

inline VECTOR operator +(const VECTOR &v)
{
    return VECTOR(v.x, v.y, v.z);
}

inline VECTOR operator -(const VECTOR &v)
{
    return VECTOR(-v.x, -v.y, -v.z);
}

inline VECTOR operator +(const VECTOR &v1, const VECTOR &v2)
{
    return VECTOR(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z);
}

inline VECTOR operator -(const VECTOR &v1, const VECTOR &v2)
{
    return VECTOR(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z);
}

inline VECTOR operator *(float f, const VECTOR &v)
{
    return VECTOR(f*v.x, f*v.y, f*v.z);
}

inline VECTOR operator *(const VECTOR &v, float f)
{
    return VECTOR(f*v.x, f*v.y, f*v.z);
}

inline VECTOR operator *(const VECTOR &v1, const VECTOR &v2)
{
    return VECTOR(v1.x*v2.x, v1.y*v2.y, v1.z*v2.z);
}

inline VECTOR operator /(const VECTOR &v, float f)
{
    float oneOverF = 1.0f / f;
    return VECTOR(oneOverF*v.x, oneOverF*v.y, oneOverF*v.z);
}

inline VECTOR operator /(const VECTOR &v1, const VECTOR &v2)
{
    return VECTOR(v1.x/v2.x, v1.y/v2.y, v1.z/v2.z);
}

inline VECTOR::operator +=(const VECTOR &v)
{
    *this = *this + v;
}

inline VECTOR::operator -=(const VECTOR &v)
{
    *this = *this - v;
}

inline VECTOR::operator *=(float f)
{
    *this = *this * f;
}

inline VECTOR::operator /=(float f)
{
    *this = *this / f;
}

inline float SquareMagnitude(const VECTOR &v)
{
    return v.x*v.x + v.y*v.y + v.z*v.z;
}

inline float Magnitude(const VECTOR &v)
{
    return sqrtf(SquareMagnitude(v));
}

inline VECTOR Normalize(const VECTOR &v)
{
    float mag = Magnitude(v);
    return v / mag;
}

inline float DotProduct(const VECTOR &v1, const VECTOR &v2)
{
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
}

inline VECTOR CrossProduct(const VECTOR &v1, const VECTOR &v2)
{
    VECTOR v;
    v.x = v1.y*v2.z - v1.z*v2.y;
    v.y = v1.z*v2.x - v1.x*v2.z;
    v.z = v1.x*v2.y - v1.y*v2.x;
    return v;
}

inline VECTOR SphericalVector(float theta, float phi)
{
    float ctheta, stheta, cphi, sphi;

    ctheta = cosf(theta);
    stheta = sinf(theta);
    cphi = cosf(phi);
    sphi = sinf(phi);

    return VECTOR(cphi * ctheta, cphi * stheta, sphi);
}

inline void GetSphereAngles(const VECTOR &v, float *theta, float *phi)
{
    *theta = atan2f(v.y, v.x);
    *phi = atanf(v.z / sqrtf(v.x*v.x + v.y*v.y));
}

inline void GetSphereAngles_Normalized(const VECTOR &v, float *theta, float *phi)
{
    *theta = atan2f(v.y, v.x);
    *phi = asinf(v.z);
}


struct HVECTOR {
    float x, y, z, w;

    HVECTOR() {}
    HVECTOR(float _x, float _y, float _z, float _w = 1.0f) {
        x = _x; y = _y; z = _z; w = _w;
    }
    HVECTOR(VECTOR v, float _w = 1.0f) {
        x = v.x; y = v.y; z = v.z; w = _w;
    }
};

inline HVECTOR operator +(const HVECTOR &v1, const HVECTOR &v2)
{
    return HVECTOR(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z, v1.w+v2.w);
}

inline HVECTOR operator -(const HVECTOR &v1, const HVECTOR &v2)
{
    return HVECTOR(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z, v1.w-v2.w);
}

inline HVECTOR operator *(float f, const HVECTOR &hv)
{
    return HVECTOR(f*hv.x, f*hv.y, f*hv.z, f*hv.w);
}

inline HVECTOR operator *(const HVECTOR &hv, float f)
{
    return HVECTOR(f*hv.x, f*hv.y, f*hv.z, f*hv.w);
}

inline VECTOR Reduce(const HVECTOR &hv)
{
    return VECTOR(hv.x, hv.y, hv.z);
}

inline VECTOR Project(const HVECTOR &hv)
{
    return VECTOR(hv.x, hv.y, hv.z) / hv.w;
}

inline float DotProduct4(const HVECTOR &v1, const HVECTOR &v2)
{
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w;
}


struct MATRIX {
    float m[4][4]; // [column][row]

    MATRIX() {}
};

inline MATRIX Identity(void)
{
    MATRIX result;

    result.m[0][0] = 1.0f;
    result.m[0][1] = 0.0f;
    result.m[0][2] = 0.0f;
    result.m[0][3] = 0.0f;
    result.m[1][0] = 0.0f;
    result.m[1][1] = 1.0f;
    result.m[1][2] = 0.0f;
    result.m[1][3] = 0.0f;
    result.m[2][0] = 0.0f;
    result.m[2][1] = 0.0f;
    result.m[2][2] = 1.0f;
    result.m[2][3] = 0.0f;
    result.m[3][0] = 0.0f;
    result.m[3][1] = 0.0f;
    result.m[3][2] = 0.0f;
    result.m[3][3] = 1.0f;

    return result;
}

inline MATRIX RotateAxis(float theta, VECTOR axis)
{
    MATRIX result;
    float s, c, ab, bc, ca, tx, ty, tz;

    s = sinf(theta);
    c = cosf(theta);
    ab = axis.x * axis.y * (1.0f - c);
    bc = axis.y * axis.z * (1.0f - c);
    ca = axis.z * axis.x * (1.0f - c);
    tx = axis.x * axis.x;
    ty = axis.y * axis.y;
    tz = axis.z * axis.z;

    result.m[0][0] = tx + c * (1.0f - tx);
    result.m[0][1] = ab + axis.z * s;
    result.m[0][2] = ca - axis.y * s;
    result.m[0][3] = 0.0f;
    result.m[1][0] = ab - axis.z * s;
    result.m[1][1] = ty + c * (1.0f - ty);
    result.m[1][2] = bc + axis.x * s;
    result.m[1][3] = 0.0f;
    result.m[2][0] = ca + axis.y * s;
    result.m[2][1] = bc - axis.x * s;
    result.m[2][2] = tz + c * (1.0f - tz);
    result.m[2][3] = 0.0f;
    result.m[3][0] = 0.0f;
    result.m[3][1] = 0.0f;
    result.m[3][2] = 0.0f;
    result.m[3][3] = 1.0f;

    return result;
}

inline MATRIX Translate(VECTOR v)
{
    MATRIX result;

    result.m[0][0] = 1.0f;
    result.m[0][1] = 0.0f;
    result.m[0][2] = 0.0f;
    result.m[0][3] = 0.0f;
    result.m[1][0] = 0.0f;
    result.m[1][1] = 1.0f;
    result.m[1][2] = 0.0f;
    result.m[1][3] = 0.0f;
    result.m[2][0] = 0.0f;
    result.m[2][1] = 0.0f;
    result.m[2][2] = 1.0f;
    result.m[2][3] = 0.0f;
    result.m[3][0] = v.x;
    result.m[3][1] = v.y;
    result.m[3][2] = v.z;
    result.m[3][3] = 1.0f;

    return result;
}

inline MATRIX Transpose(const MATRIX &m)
{
    MATRIX result;
    result.m[0][0] = m.m[0][0];
    result.m[0][1] = m.m[1][0];
    result.m[0][2] = m.m[2][0];
    result.m[0][3] = m.m[3][0];
    result.m[1][0] = m.m[0][1];
    result.m[1][1] = m.m[1][1];
    result.m[1][2] = m.m[2][1];
    result.m[1][3] = m.m[3][1];
    result.m[2][0] = m.m[0][2];
    result.m[2][1] = m.m[1][2];
    result.m[2][2] = m.m[2][2];
    result.m[2][3] = m.m[3][2];
    result.m[3][0] = m.m[0][3];
    result.m[3][1] = m.m[1][3];
    result.m[3][2] = m.m[2][3];
    result.m[3][3] = m.m[3][3];
    return result;
}

inline HVECTOR operator *(const MATRIX &m, const HVECTOR &v)
{
    HVECTOR result;
    result.x = m.m[0][0]*v.x + m.m[1][0]*v.y + m.m[2][0]*v.z + m.m[3][0]*v.w;
    result.y = m.m[0][1]*v.x + m.m[1][1]*v.y + m.m[2][1]*v.z + m.m[3][1]*v.w;
    result.z = m.m[0][2]*v.x + m.m[1][2]*v.y + m.m[2][2]*v.z + m.m[3][2]*v.w;
    result.w = m.m[0][3]*v.x + m.m[1][3]*v.y + m.m[2][3]*v.z + m.m[3][3]*v.w;
    return result;
}

inline HVECTOR operator *(const HVECTOR &v, const MATRIX &m)
{
    HVECTOR result;
    result.x = v.x*m.m[0][0] + v.y*m.m[0][1] + v.z*m.m[0][2] + v.w*m.m[0][3];
    result.y = v.x*m.m[1][0] + v.y*m.m[1][1] + v.z*m.m[1][2] + v.w*m.m[1][3];
    result.z = v.x*m.m[2][0] + v.y*m.m[2][1] + v.z*m.m[2][2] + v.w*m.m[2][3];
    result.w = v.x*m.m[3][0] + v.y*m.m[3][1] + v.z*m.m[3][2] + v.w*m.m[3][3];
    return result;
}

inline MATRIX operator *(const MATRIX &m1, const MATRIX &m2)
{
    MATRIX m;

    m.m[0][0] = m1.m[0][0]*m2.m[0][0] + m1.m[1][0]*m2.m[0][1] +
                m1.m[2][0]*m2.m[0][2] + m1.m[3][0]*m2.m[0][3];
    m.m[0][1] = m1.m[0][1]*m2.m[0][0] + m1.m[1][1]*m2.m[0][1] +
                m1.m[2][1]*m2.m[0][2] + m1.m[3][1]*m2.m[0][3];
    m.m[0][2] = m1.m[0][2]*m2.m[0][0] + m1.m[1][2]*m2.m[0][1] +
                m1.m[2][2]*m2.m[0][2] + m1.m[3][2]*m2.m[0][3];
    m.m[0][3] = m1.m[0][3]*m2.m[0][0] + m1.m[1][3]*m2.m[0][1] +
                m1.m[2][3]*m2.m[0][2] + m1.m[3][3]*m2.m[0][3];

    m.m[1][0] = m1.m[0][0]*m2.m[1][0] + m1.m[1][0]*m2.m[1][1] +
                m1.m[2][0]*m2.m[1][2] + m1.m[3][0]*m2.m[1][3];
    m.m[1][1] = m1.m[0][1]*m2.m[1][0] + m1.m[1][1]*m2.m[1][1] +
                m1.m[2][1]*m2.m[1][2] + m1.m[3][1]*m2.m[1][3];
    m.m[1][2] = m1.m[0][2]*m2.m[1][0] + m1.m[1][2]*m2.m[1][1] +
                m1.m[2][2]*m2.m[1][2] + m1.m[3][2]*m2.m[1][3];
    m.m[1][3] = m1.m[0][3]*m2.m[1][0] + m1.m[1][3]*m2.m[1][1] +
                m1.m[2][3]*m2.m[1][2] + m1.m[3][3]*m2.m[1][3];

    m.m[2][0] = m1.m[0][0]*m2.m[2][0] + m1.m[1][0]*m2.m[2][1] +
                m1.m[2][0]*m2.m[2][2] + m1.m[3][0]*m2.m[2][3];
    m.m[2][1] = m1.m[0][1]*m2.m[2][0] + m1.m[1][1]*m2.m[2][1] +
                m1.m[2][1]*m2.m[2][2] + m1.m[3][1]*m2.m[2][3];
    m.m[2][2] = m1.m[0][2]*m2.m[2][0] + m1.m[1][2]*m2.m[2][1] +
                m1.m[2][2]*m2.m[2][2] + m1.m[3][2]*m2.m[2][3];
    m.m[2][3] = m1.m[0][3]*m2.m[2][0] + m1.m[1][3]*m2.m[2][1] +
                m1.m[2][3]*m2.m[2][2] + m1.m[3][3]*m2.m[2][3];

    m.m[3][0] = m1.m[0][0]*m2.m[3][0] + m1.m[1][0]*m2.m[3][1] +
                m1.m[2][0]*m2.m[3][2] + m1.m[3][0]*m2.m[3][3];
    m.m[3][1] = m1.m[0][1]*m2.m[3][0] + m1.m[1][1]*m2.m[3][1] +
                m1.m[2][1]*m2.m[3][2] + m1.m[3][1]*m2.m[3][3];
    m.m[3][2] = m1.m[0][2]*m2.m[3][0] + m1.m[1][2]*m2.m[3][1] +
                m1.m[2][2]*m2.m[3][2] + m1.m[3][2]*m2.m[3][3];
    m.m[3][3] = m1.m[0][3]*m2.m[3][0] + m1.m[1][3]*m2.m[3][1] +
                m1.m[2][3]*m2.m[3][2] + m1.m[3][3]*m2.m[3][3];

    return m;
}


struct RAY {
    VECTOR origin;
    VECTOR dir;
};

#endif