How to calculate the intersection point between a straight and a triangle (in python)? Introduction:
I try to program a small ray tracer application in which I have planed to reflect a  ray on the surface of a triangle like the following sketches shows:

Conditions:
Given are the vectors from the triangle vertexes A,B,C and as well the two vectors P and Q from the ray. 
For solving this problem I need firstly the intersection point between the triangle and the ray (if this exists)
I have found a interesting source, which describes these scenario really good. 
CS465 Notes: Simple ray-triangle intersection http://www.cs.cornell.edu/courses/cs465/2003fa/homeworks/raytri.pdf
In this paper is the following formula to see:

Problem:
I tried to implement these formula in the following python code:
legend (python/formula):
triangle_n = n
rP = o
rD = d
p = tA -->?
def triangle_ray_intersection(tA, tB, tC, rP, rD):
triangle_xyAC = tC - tA
triangle_xyAB = tB - tA
triangle_xyBC = tC - tB

triangle_n = np.cross(triangle_xyAB, triangle_xyAC)
delta_t = np.dot((rP - tA), triangle_n)/np.dot(triangle_n,rD)
intersection_i = rD + delta_t * rP


cross_a = np.dot(np.cross(triangle_xyAB,(intersection_i-tA)),triangle_n)
cross_b = np.dot(np.cross(triangle_xyBC,(intersection_i-tB)),triangle_n)
cross_c = np.dot(np.cross(triangle_xyAC,(intersection_i-tC)),triangle_n)

But I get everytime some wrong results. I believe that I don't understand what the author from the paper means with the variable p. I thought  he means that p is just a point on the plane (and I can take a point like the given A, B or C from the triangle) but I think thats not correct. Because I get some wrong intersection point.
p.s cross_(a,b,c) are the values which should be => 0 also the condition that the intersection point is in the triangle.
Question:
Which value/variable do I have to take for the given variable p from the formula?
How i have to implement it, that I get the correct intersection point between the straight(ray) and the triangle?
How I get this value/variable from my given arguments?
Thank's you for your suggestions!
 A: Given a triangle by three vertex $\Delta\to(p_1,p_2,p_3)$ and $ p = (x,y,z)$ we have
$$
(p_1-p_2)\times(p_3-p_2) = \vec n\\
p_0 =\frac 13(p_1+p_2+p_3)
$$
are respectively the normal vector to the plane containing $\Delta$ and a point pertaining to $\Delta$. 
so we have  $\Pi\to (p-p_0)\cdot \vec n = 0$ as the plane containing $\Delta$
Now, given a line 
$$
L\to p = q_0 + \lambda\vec w
$$
we have that $\Pi \cap L = p_i$ is obtained by solving
$$
(q_0-p_0+\lambda \vec w)\cdot\vec n = 0
$$
giving
$$
p_i = q_0-\left(\frac{(q_0-p_0)\cdot\vec n}{\vec w\cdot\vec n}\right)\vec w
$$
now $p_i\in \Delta\Rightarrow (p_1-p_2)\lambda_1 +(p_3-p_2)\lambda_2 = p_i-p_2$ such that $(\lambda_1\ge 0)\cap(\lambda_2\ge 0)$ and $(p_1-p_3)\lambda_1 +(p_2-p_3)\lambda_2 = p_i-p_3$ such that $(\lambda_1\ge 0)\cap(\lambda_2\ge 0)$

import numpy as np


def dif(u,v):
    s = []
    for i in range(len(u)):
        s.append(u[i]-v[i])
    return s


def add(u,v):
    s = []
    for i in range(len(u)):
        s.append(u[i]+v[i])
    return s


def prod(u,c):
    v = []
    for i in range(len(u)):
        v.append(u[i]*c)
    return v


def isin(p,p1,p2,p3):
    u = dif(p1,p2)
    v = dif(p3,p2)
    w = dif(p,p2)
    l1 = np.dot(u,w)
    if l1 < 0: 
        return False
    l2 = np.dot(v,w)
    if l2 < 0: 
        return False
    u = dif(p1,p3)
    v = dif(p2,p3)
    w = dif(p,p3)
    l1 = np.dot(u,w)
    if l1 < 0: 
        return False
    l2 = np.dot(v,w)
    if l2 < 0: 
        return False
    return True


p1 = [10,0,0]
p2 = [0,10,0]
p3 = [0,0,10]
q0 = [10,13,15]
w  = [1,1,1]


u  = dif(p1,p2)
v  = dif(p3,p2)
n  = list(np.cross(u,v))
p0 = prod(add(add(p1,p2),p3),1/3.)
pi = add(q0,prod(w,-np.dot(dif(q0,p0),n)/np.dot(w,n)))


if isin(pi,p1,p2,p3):
    print(pi, " internal point")
else:
    print(pi, " external point")

A: The points inside the triangle can be described as the (convex) vector sum
$$\textbf A+b\,\textbf{AB}+c\,\textbf{AC}$$
where $b\ge0,c\ge0,b+c\le1$.
The points on the (half) ray are
$$\textbf P+q\,\textbf{PQ}.$$
under the condition $q\ge0$.
So to find the intersection, you can solve the vector equation
$$\textbf A+b\,\textbf{AB}+c\,\textbf{AC}=\textbf P+q\,\textbf{PQ}$$
or
$$b\,\textbf{AB}+c\,\textbf{AC}+q\,\textbf{QP}=\textbf{AP},$$
which is equivalent to a $3\times3$ linear system. Then check that the inequalities are satisfied.
