Degrees and radians are just different units for the same quantity, angular displacement. So your question is fundamentally the same issue as whether you should use, say, meters or feet to represent a distance. You just have to convert the quantity in each case to the unit that your code expects.
It's a little confusing in this case because in some cases, they've inserted the conversion factors in the formulas for you. You can recognize these conversion factors by the appearance of $\frac{180}{\pi}$ or $\frac{\pi}{180}$, which pretty much never appears in a formula except to facilitate a conversion from radians to degrees or vice-versa (respectively). For example, in
$$\Lambda = \frac{180\lambda}{\pi}$$
$\Lambda$ and $\lambda$ are actually the same quantity, but $\Lambda$ is in degrees and $\lambda$ is in radians. So when you calculate
$$\tan^{-1}\biggl[\frac{\Omega}{\sqrt{1 - \Omega^2}}\biggr]$$
if your calculator is in radian mode or you are using a function that outputs in radians, you will get an answer in radians, i.e. $\lambda$. But if your calculator is in degree mode or you are using a function that outputs in degrees, it does the conversion (the multiplication by $\frac{180}{\pi}$) for you, and you will get an answer in degrees, i.e. $\Lambda$. You don't multiply by $\frac{180}{\pi}$ again.
For simplicity, I would suggest always using trig functions that accept arguments in radians. Then when you have a value in degrees, you can convert as necessary.
In
$$\Omega = \sin\biggl(\frac{\pi\theta}{180}\biggr)\sin(\delta) + \cos\biggl(\frac{\pi\theta}{180}\biggr)\cos(\delta)\cos(\omega)$$
the formula is written with the conversion factors from degrees to radians built in. In other words, as written, it expects $\theta$ to be in degrees, but it assumes you are using radian mode (i.e. trig functions that take arguments in radians). $\theta$ needs to be converted from degrees to radians, hence the factors of $\frac{\pi}{180}$, but $\delta$ and $\omega$ are assumed to already be in radians.
In $$\delta = \frac{23.45\pi}{180}\cos\biggl(\frac{2\pi}{365}(172-J)\biggr)$$
you will notice that the factor of $\frac{\pi}{180}$ is already there to convert degrees to radians. But what it is converting is not the output of the cosine, which is a pure number (not an angle); rather, it is the $23.45^\circ$, which is the tilt of the Earth's axis relative to its orbital plane.
This answer first of all gives a simple-minded approach which has some numbers, and then a more hairy one which doesn't.
The simple approach, with numbers
So, first of all, $1367\,\mathrm{W/m^2}$ is the solar constant: it's the measured flux of energy from the Sun at the top of atmosphere (TOA), averaged over a year. So this flux is, the TOA flux for the point on the planet where the Sun is directly overhead (all other points get less). I'll call this $G_0$.
But the Earth's orbit has some eccentricity in it, so in fact sometimes this flux is a bit higher, and sometimes it's a bit lower. To first order we could model this by saying that the flux looks like
$$G = G_0\times(1 + E\cos (2\pi y))$$
where $E$ is some fudge factor based on the known eccentricity of the planet's orbit, $y$ is the time in years ($y$ is not constrained to be an integer), with $y=0$ being chosen as the point where the Earth is closest to the Sun. Observationally, $E \approx 0.03$.
Well, perhaps we want the constant in terms of day of the year, rather than year, which would be more useful. This would then look like
$$G = G_0 \times\left(1 + E\cos \left(\frac{2\pi}{365} d\right)\right)$$
Where $d$ is the day number, and $365$ is an approximate thing here: this will be OK for a while, but it will drift.
And now you have to realise that climate scientists talk to people who build spacecraft more than they talk to people like me: they work in degrees like engineers do. And $2\pi$ radians is $360$ degrees. So, finally, we get:
$$G = G_0 \times \left(1 + E\cos \left(\frac{360}{365} d\right)\right)$$
Which is your formula, and where things are working in degrees.
(I'm actually really disappointed now: I started this reply thinking 'aha, this is because you are using a model with a 360-day year (which many climate models have done historically, and many still do in fact) and I can explain this bit of obcsurity'. But no, sadly.)
A more hairy approach, without numbers
First of all we know the Sun looks quite like a black body at some temperature $T$, so the flux leaving the Sun is $\sigma T^4$ in the normal way. The total power passing through any surface surrounding the Sun is constant, so the flux at a radius $R$ is given by
$$\sigma T^4 \left(\frac{R_0}{R}\right)^2$$
Where $R_0$ is the radius of the Sun.
We could plug numbers for $T$, $R_0$ and $R$, the average radius of the Earth's orbit, into this and we will get $1367\,\mathrm{W/m^2}$. But the thing to know is how this varies with $R$, since Earth's orbit has some eccentricity. So we want to expand
$$\sigma T^4 \left(\frac{R_0}{R + \delta R}\right)^2$$ in terms of $\delta R$:
$$\begin{align}
\sigma T^4 \left(\frac{R_0}{R + \delta R}\right)^2 &= \sigma T^4\left(\frac{R_0}{R}\right)^2 \times \left(1 - 2\frac{\delta R}{R} + O(\delta R^2)\right)\\
&\approx \sigma T^4\left(\frac{R_0}{R}\right)^2 \times \left(1 - 2\frac{\delta R}{R}\right)
\end{align}$$
And now, well, we know that $\delta R$ is a periodic function of time, with the period being a year (to a good approximation & ignoring orbital variation), and there's no constant term. So, expressing time, $t$ in years, we can write $\delta R$ as
$$\delta R = \sum\limits_{n=1}^{\infty} a_n \sin(2\pi n t) + b_n \cos(2\pi n t)$$
And to first order, and adjusting the zero of $t$ suitably, we get
$$\delta R \approx a \cos(2\pi t)$$
where $a= a_1$. So plugging this into the above, we get
$$\sigma T^4\left(\frac{R_0}{R}\right)^2 \times \left(1 - 2\frac{a \cos (2\pi t)}{R}\right)$$
Whic is what we wanted.
Best Answer
The Sun will never return to the 'same spot' (e.g. relative to the Solar System centre of mass) because the orbital periods of Solar System bodies are not all rational multiples of each other. By far the largest influence on the Sun's position is caused by Jupiter (can check this by computing the force exerted by each planet). So to a decent approximation, the period of the Sun's year as you've defined it is 1 jovian year.