refract Subroutine

    subroutine refract(I, N, eta)
    !! get vector of refracted photon

        !> incident vector
        type(vector),  intent(INOUT) :: I
        !> normal vector
        type(vector),  intent(IN)    :: N
        !> \(\eta = \frac{n_1}{n_2}\)
        real(kind=wp), intent(IN)    :: eta

        type(vector)  :: T, Ntmp
        real(kind=wp) :: c1, c2

        Ntmp = N

        c1 = (Ntmp .dot. I)
        if(c1 < 0._wp)then
            c1 = -c1
        else
            Ntmp = (-1._wp) * N
        end if
        c2 = sqrt(1._wp - (eta)**2 * (1._wp-c1**2))

        T = eta*I + (eta * c1 - c2) * Ntmp 

        I = T

    end subroutine refract