Non-symmetrical DFT of real valued signal
Your choosen example frequency $\omega_0 = \frac{5 \pi} {32}$ touched a numerically sensitive case, in which the imaginary part of the theoretical DFT which should be zero is not observed so, on the numerical computation. Since the imaginary part is computed as non-zero, the phase which should be zero according to $\phi[k] = \tan^{-1}( \frac{\Im\{X_1[k]\}}{\Re\{X_1[k]\}} )$ whenever $\Im\{X_1[k]\}$ is zero, turns out to be some erratic function of $k$ (Note also that real part of $X_1[k]$ is also not computed to be exactly zero for those $k$ for which it should be, which also contributes to the erratic phase computation behaviour)

Consider the discrete time signal $x_1[n] = \cos(\frac{5\pi n}{32})$ and its 64-point DFT $X_1[k]$ (for ineteger $k \in [0:63]$) which is computed as: 

$$ X_1[k] = \sum_{n=0}^{63} \cos(\frac{5\pi n}{32}) e^{-j \frac{2\pi}{64} nk} $$ 

expand the cosine into complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{5\pi}{32} n} + e^{-j \frac{5\pi}{32} n}}) e^{-j \frac{2\pi}{64} nk} $$ 

adjust the cosine based complex exponential arguments: 
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{10\pi}{64} n} + e^{-j \frac{10\pi}{64} n}}) e^{-j \frac{2\pi}{64} nk} $$ 

separate the sums:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} {e^{j \frac{10\pi}{64} n} e^{-j \frac{2\pi}{64} nk} + 0.5 \sum_{n=0}^{63}e^{-j \frac{10\pi}{64} n}} e^{-j \frac{2\pi}{64} nk} $$ 

Now replacing numerical value of $N = 64$ with $N$ we shall see the summations:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{10\pi}{N} n} e^{-j \frac{2\pi}{N} nk} + 0.5 \sum_{n=0}^{N-1} e^{-j \frac{10\pi}{N} n}} e^{-j \frac{2\pi}{N} nk} $$

Merge the complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{2\pi n (5 -k)}{N}}+ 0.5 \sum_{n=0}^{N-1}e^{-j \frac{2\pi n (5+k)}{N}}} $$

Now the above sums add up to $N/2$ whenever $k=5$ or $k=-5$ ( this is also for $k = -5+ 64 = 59$) and they become zero for any other integer $k$ in the range $k \in [0:63]$

As you can now see, $X_1[k]$ is a real and even symmetric function of k. And its imaginary part is zero, however in your loop computation this
last condition is not realized due to numerical issues and therefore you observe erratic phase computation behaviour.