fresnel Function

    function fresnel(I, N, n1, n2) result (tir)
    !! calculates the fresnel coefficents

        use ieee_arithmetic, only : ieee_is_nan

        !> reffractive indicies
        real(kind=wp), intent(IN) :: n1, n2
        !> incident vector
        type(vector),  intent(IN) :: I
        !> Normal vector
        type(vector),  intent(IN) :: N

        real(kind=wp) :: costt, sintt, sint2, cost2, tir, f1, f2

        costt = abs(I .dot. N)

        sintt = sqrt(1._wp - costt * costt)
        sint2 = n1/n2 * sintt
        if(sint2 > 1._wp)then
            tir = 1.0_wp
            return
        elseif(costt == 1._wp)then
            tir = 0._wp
            return
        else
            sint2 = (n1/n2)*sintt
            cost2 = sqrt(1._wp - sint2 * sint2)
            f1 = abs((n1*costt - n2*cost2) / (n1*costt + n2*cost2))**2
            f2 = abs((n1*cost2 - n2*costt) / (n1*cost2 + n2*costt))**2

            tir = 0.5_wp * (f1 + f2)
        if(ieee_is_nan(tir) .or. tir > 1._wp .or. tir < 0._wp)print*,'TIR: ', tir, f1, f2, costt,sintt,cost2,sint2
            return
        end if
    end function fresnel