struct Ellipsoid : public Intersectable {
    vec3 center;
    float a;
    float b;
    float c;

    Ellipsoid(const vec3& _center, float _a, float _b, float _c, Material* _material) {
        center = _center;
        a = _a;
        b = _b;
        c = _c;
        material = _material;
    }

    Hit intersect(const Ray& ray) {
        Hit hit;
        mat4 M(
            1.0f / (a*a), 0, 0, 0,
            0, 1.0f / (b*b), 0, 0,
            0, 0, 1.0f / (c*c), 0,
            0, 0, 0, -1
        );

        vec4 v(ray.dir, 0);
        vec4 eye(ray.start-center, 1);

        float A = dot(v*M,v);

        float B = dot(eye*M, v) + dot(v*M, eye);

        float C = dot(eye*M,eye);

        float discr = B * B - 4.0 * A * C;
        if (discr < 0) return hit;
        float sqrt_discr = sqrtf(discr);
        float t1 = (-B + sqrt_discr) / 2.0 / A;
        float t2 = (-B - sqrt_discr) / 2.0 / A;
        if (t1 <= 0 && t2 <= 0) return hit;
        if (t1 <= 0 && t2 > 0)       hit.t = t2;
        else if (t2 <= 0 && t1 > 0)  hit.t = t1;
        else if (t1 < t2)            hit.t = t1;
        else                         hit.t = t2;

        hit.position = ray.start + ray.dir * hit.t;
        float x_c = hit.position.x - center.x;
        float y_c = hit.position.y - center.y;
        float z_c = hit.position.z - center.z;
        hit.normal = vec3(x_c/(a*a), y_c/(b*b), z_c/(c*c)).normalize();
        if (dot(hit.normal, ray.dir) > 0) hit.normal = hit.normal * (-1);
        hit.material = material;
        return hit;
    }
};