I am building a simplified model solar system in GeoGebra. Celestial objects are placed in a heliocentric coordinate system with the sun at the origin, the x-y-plane as the ecliptic, and the x-axis aligned with Earth's March Equinox. Earth's center is located at the point (EarthOrbitRadius cos(2π t), EarthOrbitRadius sin(2π t)), where t is the time in years. Earth's equator and axis of rotation are rotated by AxialTilt° about the line y = EarthOrbitRadius sin(2π t), which marks the directions of the March and September equinoxes.
I want to know the right ascension (RA, the angle eastward along the celestial equator between the line marking the March equinox and the hour circle of the target point, in the range 0° ≤ RA < 360°) and declination (DEC, the angle along the hour circle of the target point, perpendicular to the celestial equator, in the range −90° ≤ DEC ≤ 90° with negative values being south of the celestial equator and positive values being north) for the sun as measured from Earth. The hour circle is the great circle that passes through the target point and the two celestial poles; in the case of the Earth, these are the north and south poles along the axis of rotation. Put another way, the hour circle of a point lies on the plane formed between the line of the axis of rotation and the line through the point and the center of the Earth.
If the Earth's axis was not tilted, this would be easy. The RA of the sun would be simply 2π t, and the DEC of the sun would be a constant 0°. With the axis tilted, however, the values change. The RA acquires a wobble, and the DEC becomes approximately (but not exactly, it also has a wobble) AxialTilt sin(2π t).
The wobbles in the RA and DEC values appear to be roughly sinusoidal, but I can't figure out an exact form for them. The RA anomaly looks approximately like sin(4π t) with the peaks shifted and scaled down vertically by some relation to the AxialTilt value, while the DEC anomaly looks approximately like −sin(2π t) − sin(6π t), also scaled vertically by some relation to the AxialTilt value. In both cases, the larger the AxialTilt value, the larger the wobble.
Is there an exact form for these wobbles?
The solution to this problem can be found relatively easily by transforming to the coordinate system of the equatorial plane in a geocentric manner. We know that the equatorial plane is generated by the normal vector to it, which is the axis of rotation by definition. We need to find two vectors to span the plane. The equinox positions offer us a good place to start, because there the normal vector $\hat{X}$ pointing from the sun towards the equator belongs to both the orbital and equatorial plane. We assume thus that in our geocentric system the Sun rotates around the Earth and that at time $t=0$ it is to be found in the autumn equinox. In the orbital coordinate system (O) the sun has coordinates
$$r_s=R\cos2\pi t ~\hat{x}+R\sin 2\pi t ~\hat{y}$$
and also $\hat{X}=\hat{x}$. We know the approximately constant tilt angle $\alpha$ of the Earth's axis so we can assume that
$$\hat{n}=(\sin\alpha\cos\beta, \sin\alpha\sin\beta, \cos\alpha)$$
for some angle $\beta$ that parametrizes where the axis is pointing at the equinox. We want to use this vector as part of our equatorial orthornormal basis. Demanding that $\hat{n}\cdot\hat{X}=0$ requires that $\beta=\pi/2 , 3\pi/2$, and we pick $\beta=\pi/2$ without loss of generality. Taking the vector $\hat{Y}=\hat{n}\times \hat{X}$ as the third basis vector we finally see that the basis vectors can be written in (O) coordinates:
$$\hat{X}=(1,0,0),\hat{Y}=(0,\cos\alpha, -\sin\alpha),\hat{n}=(0, \sin\alpha, \cos\alpha)$$
In equatorial coordinates we can compute the Sun's position, which reads
$$\frac{r_s}{||r_s||}=(\cos 2\pi t, \sin 2\pi t\cos \alpha, \sin 2\pi t \sin\alpha)$$
We readily read off the declination and right ascension angles from the formula
$$\frac{r_s}{||r_s||}=(\cos(DEC)\cos(RA),\cos(DEC)\sin(RA), \sin(DEC))$$
by comparing to the vector we found:
$$DEC=\arcsin\left(\sin\alpha \sin 2\pi t \right)~, ~t\in(0,1)~~~~ RA= \begin{Bmatrix}\arctan\left(\cos\alpha \tan 2\pi t \right)~,~t\in (0,1/2)\\ \pi+\arctan\left(\cos\alpha \tan 2\pi t)~,~ t\in (1/2,1 \right)\end{Bmatrix}$$