Does anyone know a general way to find all the vectors at a given angle a from a vector?
In two dimensions, it is one of these vectors (x,y) with parameter r. $$(- y \sin(a) + x \cos(a), x \sin(a) + y \cos(a)) r$$ $$(y \sin(a) + x \cos(a), - x \sin(a) + y \cos(a)) r$$
What is the formula for three dimensions? In this case, the vectors would all lie in a cone around the main vector, so there should be two parameters. Would spherical coordinates help?
Is there a formula for n dimensions?
For 3 and higher dimensions, it is the same.
For any vector $y:=(y_1,\dots,y_n)$ that has certain angle with the given vector $x:=(x_1,\dots,x_n)$, obviously $x,y$ span a plane. You can always use an orthogonal transformation, which is angle preserving, to bring the two into the plane $(x_1,x_2,0,\dots,0)$. So you are back in the old game!
In this way, it is not difficult to see that all such $y$'s are generated from a particular one by orthogonal transformations fixing your original vector $x$. The parameter space looks like $SO(n-1)\times\mathbb R$, i.e. direction and length.
The actual carry out is relatively easy. Just complete $x$ to a complete basis. Perform Gram-Schmidt orthogonalization to get an orthonormal basis $x_1,\dots,x_n$, with $x$ parallel to $x_1=(1,0,\dots,0)$. Assume w.l.o.g. that $x=x_1$ Now construct your $y=x_1\cos\theta+x_2\sin\theta=(\cos\theta,\sin\theta,0,\dots,0)$ under the new basis, where $\theta$ is the desired angle $y$ makes with $x$. Elements of $SO(n-1)$, i.e. rotation matrix preserving $x_1$ should have the following form: $$ A=\left[\begin{array}{c|c} 1 & 0\\ \hline 0 & A' \end{array}\right] $$ where $A'$ is an $(n-1)\times (n-1)$ special (with determinant $1$) orthogonal matrix. This is the matrix you need to rotate $y$ to get other admissible vectors.