Combining rotations for orbital and rotational motion of a planet If I have a planet orbiting the Sun (assuming circular orbit) at angular velocity $\Omega$ and rotating about its axis at $\omega$. I also have a normal to the surface of the planet $\vec{n}_{\rm surf}$, the normal to the orbital plane $\vec{n}_{\rm o}$, and the direction to the Sun from the planet and a normal to the planet's surface $\vec{n}_{\rm sun}$. The question is how can I get the dot product $\vec{n}_{\rm sun}(t)\cdot\vec{n}_{\rm surf}(t)$ at any time during the orbit given these values at some $t=0$, eg. in perihelion? The problem is that these are two rotations that have to be combined. I have two approaches:
My first approach is to be in the planet's system of coordinates, see the picture:

The problem here is that $\vec{n}_{\rm o}$ is changing due to the rotation about $z$-axis of the planet. So I am thinking:


*

*Rotate $\vec{n}_{\rm o}(0)$ about $\vec{z}$ by ${\rm d}t * \omega$ to get $\vec{n}_{\rm o}(1)$ 

*Rotate $\vec{n}_{\rm sun}(0)$ about $\vec{n}_{\rm o}(0)$ by ${\rm d}t * \Omega$ to get $\vec{n}_{\rm sun}(1)$

*Rotate $\vec{n}_{\rm o}(1)$ about $\vec{z}$ by ${\rm d}t * \omega$ to get $\vec{n}_{\rm o}(2)$ 

*Rotate $\vec{n}_{\rm sun}(1)$ about $\vec{n}_{\rm o}(1)$ by ${\rm d}t * \Omega$ to get $\vec{n}_{\rm sun}(2)$

*Continue, until I get $\vec{n}_{\rm sun}(t)$

*Because $\vec{z}$ is constant, I can get $\vec{n}_{\rm surf}(t)$ directly as rotation of $\vec{n}_{\rm surf}(0)$ about $z$-axis by $t * \omega$


I know this is only an approximate solution and quality depends on ${\rm d}t$ step. So I was thinking, would I not get the same results by this second approach.
I assume that I am in the coordinates of the sun (neglecting its rotation), such that the $z$-axis is normal to the orbital plane, $x$ axis is such that it points to the perihelion. See the picture:

Now my approach would be:


*

*Rotate $\vec{n}_{\rm sun}(0)$ about $z$ axis by $t*\Omega$ to get $\vec{n}_{\rm sun}(t)$.

*Rotate $\vec{n}_{\rm surf}(0)$ about the $\vec{n}_{\rm spin}$ by $t*\omega$ to get $\vec{n}_{\rm surf}(t)$.


The reason why I think I can do this is that $\vec{n}_{\rm spin}$ will not change by orbital rotation and neither will $\vec{n}_{\rm surf}$. But I am not really sure about this assumption.
 A: The dot product $\vec{n}_{sun}(t) \cdot \vec{n}_{surf}(t)$ can be calculated as $||\vec{n}_{sun}|| \cdot ||\vec{n}_{surf}|| \cos(\theta(t))$.  Note that I carried the time variable only to the angle as I am assuming from your description that the magnitude of the normal vectors is constant in time.  From that point use your angles $\Theta = \Theta_o + \Omega \cdot t$ and $\theta = \theta_o + \omega \cdot t$ and combine this to construct the Euler angles (https://en.wikipedia.org/wiki/Euler_angles) of the planet.
This should get you on your way.  For the details I would put on my professor hat and leave it to the interested reader.
A: (I'm not quite sure that I understand the assumptions of your problem correctly.)
"The reason why I think I can do this is that $\vec{n}_{spin}$ will not change by orbital rotation and neither will $\vec{n}_{surf}$. But I am not really sure about this assumption" 
I think this assumption is somewhat lying in the definition of $\vec{\Omega}$ and $\vec{\omega}$. [$\vec{\Omega}$ just moves the planet solidly around the orbit (and its value determines the angular velocity of the planet's center around the sun), and $\vec{\omega}$ determines the value of the angular velocity of the planet's surface around an axis (call it $z'$), which is fixed as seen by an observer living on the planet.]  
As you would probably know, if a vector with fixed length rotates around an axis with an angular velocity, the cross product of the angular velocity and that vector at any specific time, gives the rate of change of that vector at that time. So we have: 
$$ \frac{d \vec{n}_{sun}}{d t} = \vec{\Omega}\times\vec{n}_{sun},$$ 
and also:
$$\frac{d \vec{n}_{surf}}{d t} = \vec{\omega} \times \vec{n}_{surf},$$
then we can proceed to the inner product:
\begin{align}
 \frac{d}{dt} (\vec{n}_{sun}. \vec{n}_{surf}) &= (\vec{\Omega}\times \vec{n}_{sun}).\vec{n}_{surf} + \vec{n}_{sun}.(\vec{\omega} \times \vec{n}_{surf})\\ &= \vec{n}_{surf}.[(\vec{\Omega}-\vec{\omega})\times \vec{n}_{sun}]
\end{align}
And then choose a suitable coordinate and solve this equation (Actually it seems that solving for $\vec{n}_{surf}$ and  $\vec{n}_{sun}$ straightly is way easier)
Hope that it helps.
A: You have four different coordinate systems at play here:


*

*Coordinates fixed to earth. Think latitude/longitude/height, or some cartesian system where London stays at a single constant coordinate.

*Coordinates fixed to the earth's center of mass, but rotationally fixed to the celestial sphere. In this coordinate system, places on the equator move eastwards at a speed of $\frac{40000km}{24h} = 1667\frac{km}{h}$.

*Coordinates fixed to the sun's center of mass, but rotationally fixed to the celestial sphere. In this system, the earth's poles are not stationary anymore, they move around the sun at $\frac{150\cdot10^6km\cdot2\pi}{365.25\cdot24h} \approx 100000km/h$

*Coordinates fixed to the sun's center of mass, but rotationally fixed to follow earth's center of mass on its yearly journey around the sun. In this system, the earth's axis itself rotates as the year goes by.
Between these three coordinate systems, you have simple transformations:


*

*1 <-> 2: Rotate about earth's axis.
Both, normal and location vectors are translated.

*2 <-> 3: Translate along the vector connecting earth's and the sun's center of mass.
This transformations does not change any normal vectors.

*3 <-> 4: Rotate about the sun within earth's orbital plane.
Again, both normals and location vectors are transformed.
All these transformations are time dependent.

Your vector $\vec{n}_{surf}$ is a constant in the coordinate system 1; $\vec{n}_{sun}$ is a constant in system 4.
As such, all you need to do is, either transform $vec{n}_{surf}$ until it is expressed in the coordinates of system 4, or transform $vec{n}_{sun}$ until it is expressed in the coordinates of system 1. Either way, once you have both vectors in the same coordinate system, calculating their dot product is trivial.
Notes:


*

*The transformations are time dependent, so after the first transformation the vectors will be functions of time $\vec{n}_{surf}(t)$ and  $\vec{n}_{sun}(t)$.

*The second transformation is a noop on normal vectors. All that matters are the two rotations. You simply need to apply two time dependent rotations around two different axes.
A: I assume that the planet is  rotating about the surface axis $\vec{n}_f$ with the rotation angle $\omega\,t$, and  rotate about the sun with the rotation axis $\vec{n}_0$ with the angle $\Omega\,t$.
Where $\vec{n}_f$ is your $\vec{n}_{surf}$ vector
Thus:
$$R_p=R_{n0}(\Omega\,t,\vec{\hat{n}}_0)\,R_f(\omega\,t,\vec{\hat{n}}_f)$$
Where $R_{n0}\,,R_p$ are Rodriguez' Rotation matrices and $R_p$ is the rotation matrix between the planet fixed frame and inertial frame.
The vector $\vec{n}_{sun}$ is rotating with $R_{n0}$ and the vector $\vec{n}_f$ is rotating with $R_p$ matrix, so taking the dot product:
with :
$$\vec{n}_s=R_{n0}\,\vec{n}_s(0)$$
$$\vec{n}_f=R_{p}\,\vec{n}_f(0)$$
Thus:
$$\vec{n}_{s}\cdot \vec{n}_{f}=(\vec{n}_s(0))^T\,(R_{n0})^T\,R_{n0}\,R_f\,\vec{n}_f(0)
=(\vec{n}_s(0))\cdot \left[\,R_f(\omega\,t,\vec{\hat{n}}_f)\,\vec{n}_f(0)\right]$$
