Trying to construct an elliptic curve pairing (Tate) example. Can someone help me figure out where my mistake is? I'm trying to construct an example of the Tate pairing, but for some reason I'm not getting the bilinearity property that I should be getting. Here's my setup:
Take $q=29$ and $$E/\mathbb{F}_q : y^2 = x^3 + 21x $$
$E$ has $26$ points, and using $r = 13$ for the $r$-torsion, we have embedding degree $k = 3$.
Pick a point in the $r$-torsion, $P = (28,6)$.
Next we extend $E$ to $E/\mathbb{F}_{q^k}$ by creating 
$$\mathbb{F}_{q^k} := \mathbb{F}_q(u)$$
where $$u^3+u+4 \equiv 0$$
because that polynomial is irreducible mod $29$.
$E/\mathbb{F}_{q^k}$ has $24674$ elements, so take $h=\frac{24674}{r^2} = 146$. Next we pick any random point on $E/\mathbb{F}_{q^k}$, call it $Q$. We can send $Q$ to the trace-zero subgroup of the $r$-torsion by multiplying it by $h$ (this brings it into the $r$-torsion), and then applying the anti-trace map
$$Q := k*Q - \text{Trace}(Q)$$

From here, I'm just using a code for Miller's algorithm that I found here. The code's example works. My issue is that when I compute the Tate pairing for $P$ and $Q$, I don't get bilinearity...
I'm happy to post my full code if anyone's willing to take a look. It's written in Magma.

Alex, there is a problem in all the magma codes posted in Pairings by Craig Costello. 
I bet you didn't get the bilinearity because in fADD function there is a mistake. Change this line v:=F!(x-(lambda^2-P[1]-P[2])); with v:=F!(x-(lambda^2-P[1]-Q[1]));
That's the correct way to get the vertical line equation. Unfortunately the same mistake is present in all the other magma code samples in that paper.
Btw, I've spent 2 days to figure it out that mistake .. But, I am happy now.
