intersectEllipse Function

    logical function intersectEllipse(orig, dir, t, centre, semia, semib)
    !! calculates where a line, with origin:orig and direction:dir hits a ellipse, centre:centre and axii:semia, semib
    !! returns true if intersection exists
    !! returns t, the paramertised parameter of the line equation
    !! adapted from scratchapixel and pbrt
    !! need to check z height after moving ray
    !! if not this is an infinte ellipse-cylinder
    !! ellipse lies length ways along z-axis
    !! semia and semib are the semimajor axis which are the half width and height.

        !> Direction vector of the ray
        type(vector),  intent(IN)  :: dir
        !> origin of the ray
        type(vector),  intent(IN)  :: orig
        !> Centre of the ellipse
        type(vector),  intent(IN)  :: centre
        !> distance from orig to the intersection point
        real(kind=wp), intent(OUT) :: t
        !> Half width of the ellipse
        real(kind=wp), intent(IN)  :: semia
        !> Half height of the ellipse
        real(kind=wp), intent(IN)  :: semib

        type(vector)  :: L
        real(kind=wp) :: t0, t1, a, b, c, tmp, semia2div, semib2div

        intersectEllipse = .false.

        semia2div = 1._wp / semia**2
        semib2div = 1._wp / semib**2

        L = orig - centre
        a = semia2div * dir%z**2 + semib2div * dir%y**2
        b = 2._wp * (semia2div * dir%z * L%z + semib2div * dir%y * L%y)
        c = semia2div * L%z**2 + semib2div * L%y**2 - 1._wp

        if(.not. solveQuadratic(a, b, c, t0, t1))return
        if(t0 > t1)then
            tmp = t1
            t1 = t0
            t0 = tmp
        end if
        if(t0 < 0._wp)then
            t0 = t1
            if(t0 < 0._wp)return
        end if

        t = t0
        intersectEllipse = .true.
        return
    end function intersectEllipse