Matrix3 ComputeJacobiRotation(const Matrix3& matrix)
{
    std::vector<int> pq;

    if (abs(matrix[0][1]) > abs(matrix[0][2]))
    {
        if (abs(matrix[0][1]) > abs(matrix[1][2]))
            pq = { 0, 1 };
        else
            pq = { 1, 2 };
    }
    else
    {
        if (abs(matrix[0][2]) > abs(matrix[1][2]))
            pq = { 0, 2 };
        else
            pq = { 1, 2 };
    }

    // [c, s]
    // [-s, c]

    float beta = (matrix[pq[1]][pq[1]] - matrix[pq[0]][pq[0]]) /
                 (2 * matrix[pq[0]][pq[1]]);

    float tan = (beta > 0.0f ? 1.0f : -1.0f) / (abs(beta) + sqrt((beta * beta) + 1.0f));

    float cos2 = 1.0f / ((tan * tan) + 1);
    float sin2 = 1.0f - cos2;

    float cos = sqrt(cos2);
    float sin = sqrt(sin2);

    Matrix3 result;

    result.SetIdentity( );

    result[pq[0]][pq[0]] = cos;
    result[pq[0]][pq[1]] = sin;
    result[pq[1]][pq[1]] = cos;
    result[pq[1]][pq[0]] = -sin;

    return result;
}