#ifndef HEADER_RAY_H__
#define HEADER_RAY_H__

#include "v.h"

template<class T> class Ray
{
private:
        bool f_intersected;
public:
        V<T> r;
        V<T> d;

        Ray(const T& x1, const T& y1, const T& z1, const T& x2, const T& y2, const T& z2) : r(x1, y1, z1), d(x2, y2, z2)
        {
                f_intersected = true;
        }

        bool isIntersected()
        {
                return f_intersected;
        }

        V<T> xyPlaneIntersection(const T& th)
        {
                V<T> v;
                T dist = th - r.z;
                v.x = r.x + dist * d.x / d.z;
                v.y = r.y + dist * d.y / d.z;
                v.z = th;
                return v;
        }

        bool intersectsRectangle(const T& th, const T& a, const T& b)
        {
                V<T> v = xyPlaneIntersection(th);
                bool intrsct = (v.x >= -a / 2.0) && (v.x <= a / 2.0) && (v.y >= -b / 2.0) && (v.y <= b / 2.0);
                f_intersected = f_intersected && intrsct;
                return intrsct;
        }

        bool intersectsCircle(const T& th, const T&  r)
        {
                V<T> v = xyPlaneIntersection(th);
                auto tr = v.x * v.x + v.y * v.y;
                bool intrsct = tr <= r * r;
                f_intersected = f_intersected && intrsct;
                return intrsct;
        }

        template<class Ang> bool intersectsSphere(const T& tH, const T& tR, const Ang& salf, const Ang& sbet)
        {
                T l = tH + tR;
                T fx = l * sin(salf) / sqrt(1.0 + tan(sbet)*tan(sbet) * cos(salf)*cos(salf));
                T fy = l * sin(sbet) / sqrt(1.0 + tan(salf)*tan(salf) * cos(sbet)*cos(sbet));
                double sgn = 1;
                if (cos(salf) < 0.0 || cos(sbet) < 0.0)
                {
                        sgn = -1.0;
                }
                T fz = sgn * l / sqrt(1.0 + tan(salf) * tan(salf) + tan(sbet) * tan(sbet));

                auto as = d.x * d.x + d.y * d.y + d.z * d.z;
                auto bs = 2.0 * (d.x * (r.x - fx) + d.y * (r.y - fy) + d.z * (r.z - fz));
                auto cs = (r.x - fx) * (r.x - fx) +
                              (r.y - fy) * (r.y - fy) +
                                  (r.z - fz) * (r.z - fz)- tR * tR;

                bool intrsct = false;
                auto Ds = (bs * bs - 4.0 * as * cs);
                decltype(Ds) zero(0);
                if (Ds >= zero)
                {
                        auto tau = (-bs + sqrt(Ds)) / 2.0 / as;
                        T x = r.x + tau * d.x;
                        T y = r.y + tau * d.y;
                        T z = r.z + tau * d.z;
                        intrsct = z >= T(0);
                }
                f_intersected = f_intersected && intrsct;
                return intrsct;
        }
};