You really don't need to go that painful way.
It pays to collect your $\theta_a \equiv \theta~ n_a$, where $n_a$ are the components of a unit 3-vector $\hat{\boldsymbol{n}}$.
As a consequence, for spin 1, that is, the triplet representation of SU(2), the celebrated Rodrigues formula gives
$$
e^{i\theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}=I_{3}+i(\boldsymbol{\hat
{n}}\cdot\boldsymbol{J})\sin{\theta}+(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})^{2}(\cos\theta-1)~.
$$
The reason is that the eigenvalues of $\hat{\boldsymbol{n}}\cdot\boldsymbol{J}$ are the eigenvalues of
$J_{3}=\mathrm{diag}(1,0,-1)$, and the Cayley-Hamilton theorem (dictating this matrix satisfies its own characteristic equation) reduces any higher power of it to a linear combination of the identity and the first two powers.
It is then straightforward to obtain this expression.
(You might find the compact formula for any irrep in Curtright et al, 2014.)