Are the irreps of SO(4) necessarily real? I'm not a group theory buff, but I was confused when I used some code to generate the irreducible representatives of SO(4) and found the resulting matrix elements to be complex. Is this possible, or must the irreps of SO(4) be real? (It would seem to me that the distinction of "orthogonal" ought to be reserved for real matrices, which is why this seems strange to me, but I'm probably missing something).
 A: This answer assumes comfort with representation theory, in particular knowledge of the representation theory of $SU(2)$. Try Peter Woit's notes for a good introduction to the latter.
All the irreps of $SO(4)$ can be defined over the reals. We have $SO(4) \cong \left( SU(2) \times SU(2)\right) / \langle (-\mathrm{Id}, -\mathrm{Id}) \rangle$. To repeat, we are quotienting by the two element subgroup whose nontrivial element is $ (-\mathrm{Id}, -\mathrm{Id})$.
Let $V$ be the standard, two dimensional complex representation of $SU(2)$. Let $V^p = \mathrm{Sym}^{p} V$. Every complex irrep of $SU(2)$ is of the form $V^p$ for some $p \geq 0$; every complex irrep of $SU(2) \times SU(2)$ is of the form $V^p \boxtimes V^q$. 
The element $- \mathrm{Id}$ acts on $V^p$ by $(-1)^p$. So the representation $V^p \boxtimes V^q$ factors through $SO(4)$ if and only if $p \equiv q \mod 2$.
Case 1: The indices $p$ and $q$ are even. I claim that $V^p$ can then be defined over $\mathbb{R}$ (and likewise $V^q$), so $V^p \boxtimes V^q$ can be defined over $p$. In other words, the irreducible complex representations of $SO(4)$ are $V^p \boxtimes V^q$, for $(p,q)$ even.
The representation $V^p$ factors through $SO(3)$ when $p$ is even (since $- \mathrm{Id}$ acts by $1$). We need to see that all the representations of $SO(3)$ can be defined over the reals. $V^2$ is just the standard $3$-dimensional representation of $SO(3)$, and so defined over the reals. We have 
$$\mathrm{Sym}^r V^2 = V^{2r} \oplus \bigoplus_{s<r} (V^{2s})^{\oplus (\mbox{some multiplicity we don't care about})}$$
The representation $\mathrm{Sym}^r V^2$ is defined over $\mathbb{R}$, so are all of its isotypic components are, so $V^{2r}$ is. 
Case 2: The indices $p$ and $q$ are odd. By the previous case, $V^{p-1} \boxtimes V^{q-1}$ is defined over $\mathbb{R}$. Also, $V^1 \boxtimes V^1$ is the standard $4$-dimensional representation of $SO(4)$, which is defined over $\mathbb{R}$. 
So $(V^{p-1} \boxtimes V^{q-1}) \otimes (V^1 \boxtimes V^1)$ is defined over $\mathbb{R}$.
We have
$$(V^{p-1} \boxtimes V^{q-1}) \otimes (V^1 \boxtimes V^1) \cong (V^{p-1} \otimes V^1) \boxtimes (V^{q-1} \otimes V^1) \cong (V^p \oplus V^{p-2}) \boxtimes (V^q \oplus V^{q-2})$$
$$\cong V^p \boxtimes V^q \oplus (\mbox{stuff not containing a $(p,q)$ factor}).$$
Again, we have recovered $V^{p} \boxtimes V^q$ as an isotypic component of something defined over $\mathbb{R}$, so it is defined over $\mathbb{R}$.
What might be wrong with your program? Well, either your program is giving you the wrong representations, or it is giving you representations that look complex but if you wrote them in a better basis, would be visibly real. We can probably take a good guess as to which it is if you record how you think a matrix of the form
$$\begin{pmatrix}
\cos \alpha & \sin \alpha & 0 & 0 \\
-\sin \alpha & \cos \alpha & 0 & 0 \\
0 & 0 & \cos \beta & \sin \beta  \\
0 & 0 & -\sin \beta & \cos \beta  \\
\end{pmatrix}$$
acts on one of your complex representations.
