Description
I'm trying to reproduce a Solar System mathematically using ellipses. The initial conditions are:
- All ellipses will have a common foci $(0,0)$, but different centers
- A planet orbiting an ellipse will have a period $T$
- The position of a planet in the ellipse will be the reflected point (from the center of the ellipse) given the angle $\theta$ at time $t$ (*)
- Every ellipse can have a different initial phase angle $\theta_0$
- Every ellipse can be rotated about the common foci $(0,0)$
To simplify the problem. We will assume conditions $4.$ and $5.$ are nonexistent. Here's a simple representation of a single ellipse:
where
- $(j,k)$ is the center of the ellipse
- $(x,y)$ is the result of the ellipse formula at angle $\theta$
- $(x',y')$ is the reflected point using the center as a mirror
- $\omega$ is the angle at point $(x',y')$
*: This is needed to simulate Kepler's second law
Problem
Given two ellipses (whose variables are denoted by subindex $_1$ and $_2$ in the following formulas), I'd like to find a formula that describes at which times $t$, $\omega_1$ and $\omega_2$ are equal. In other words, find the equation of $\omega$ in relation to time $t$. Here's a gif that shows the specified system of two orbits:

where
- The green dot is the center of the large ellipse
- The red dot is the center of the smaller ellipse
- The period of the large ellipse is $T_1 = 6$s
- The period of the smaller ellipse is $T_2 = 14$s
- The eccentricity of the large ellipse is 0.77
- The eccentricity of the smaller ellipse is 0.75
- The directrix of both ellipses is 2
So far, I've been able to figure out at which times $t$ both $\theta$ are equal:
We know how $\theta$ varies in relation to $t$, $$\theta = \frac{2\pi t}{T} + \theta_0$$
Since we have stablished that $\theta_0$ doesn't exist, $$\theta = \frac{2\pi t}{T}$$
We know that at time $t$, both $\theta_1$ and $\theta_2$ must be equal: $$\theta_2 = \theta_1 + 2\pi k\ \ ,\ k \in \mathbb{Z}$$
where $2\pi k$ represents the periodicity of the orbit. If we substitute: $$\frac{2\pi t}{T_2} = \frac{2\pi t}{T_1} + 2\pi k\ \ ,\ k \in \mathbb{Z}$$
Solve for t: $$t = \frac{k}{\frac{1}{T_2}-\frac{1}{T_1}}$$
This formula gives me at which time $t$ both ellipses will have the same $\theta$. If we substitute $T_1$ and $T_2$ with the periods of the ellipses in the example gif, we got for $k=1$, $t=10.5$s. If you look at the gif, after 10.5s, both red lines will be aligned, meaning $\theta_1$ and $\theta_2$ will be equal. Another way of finding this formula is using the polar expression of an ellipse:
We know that the angle $\theta$ will be $$\theta = \arctan{\frac{y}{x}}$$
And we know that the point $(x,y)$ in polar coordinates is represented as: $$(x,y) = (r\cos{\theta}, r\sin{\theta})$$
Where $r$ is the distance from the center to the point $(x,y)$, and this distance is given by the polar formula of an ellipse: $$r=\frac{ed}{1+e\cos{\theta}}$$ $$(x,y) = \left(\frac{ed}{1+e\cos{\theta}}\cos{\theta}, \frac{ed}{1+e\cos{\theta}}\sin{\theta}\right)$$
Going back to the formula described at step $1.$, we have $$\theta = \arctan{\frac{y}{x}} = \arctan{\frac{\frac{ed}{1+e\cos{\theta}}\sin{\theta}}{\frac{ed}{1+e\cos{\theta}}\cos{\theta}}}$$
Simplifying, we have $$\theta = \arctan{\frac{y}{x}} = \arctan{\frac{\sin{\theta}}{\cos{\theta}}} = \arctan{(tan{\theta})}$$
Again, we want both $\theta_1$ and $\theta_2$ to be equal: $$\arctan{\left(tan{\frac{2\pi t}{T_2}}\right)} = \arctan{\left(tan{\frac{2\pi t}{T_1}}\right)}$$
Which is equivalent to $$\frac{2\pi t}{T_2} = \frac{2\pi t}{T_1} + 2\pi k\ \ ,\ k \in \mathbb{Z}$$
So far so good. Now we want to find at which time $t$, both $\omega_1$ and $\omega_2$ are equal. We can use the second approach:
We know that the point (x',y') is $$(x',y')=2(j,k)-(x,y) = (2j-x,2k-y)$$
And we already know how to represent $x$ and $y$ in relation to time: $$(x,y) = \left(\frac{ed}{1+e\cos{\theta}}\cos{\theta}, \frac{ed}{1+e\cos{\theta}}\sin{\theta}\right)$$
This means that $(x',y')$ will be $$(x',y') = \left(2j-\frac{ed\cos{\theta}}{1+e\cos{\theta}}, 2k-\frac{ed\sin{\theta}}{1+e\cos{\theta}}\right)$$
Which means that the angle $\omega$ will be $$\omega = \arctan{\left(\frac{2k-\frac{ed\sin{\theta}}{1+e\cos{\theta}}}{2j-\frac{ed\cos{\theta}}{1+e\cos{\theta}}}\right)}$$
Since we want $\omega_1$ and $\omega_2$ to be equal $$\omega_2 = \omega_1 + 2\pi k$$ $$\arctan{\left(\frac{2k_2-\frac{e_2d_2\sin{\theta_2}}{1+e_2\cos{\theta_2}}}{2j_2-\frac{e_2d_2\cos{\theta_2}}{1+e_2\cos{\theta_2}}}\right)} = \arctan{\left(\frac{2k_1-\frac{e_1d_1\sin{\theta_1}}{1+e_1\cos{\theta_1}}}{2j_1-\frac{e_1d_1\cos{\theta_1}}{1+e_1\cos{\theta_1}}}\right)} + 2\pi k$$
and we know that $k_2=k_1=0$ since condition $5.$ doesn't apply: $$\arctan{\left(\frac{-\frac{e_2d_2\sin{\theta_2}}{1+e_2\cos{\theta_2}}}{2j_2-\frac{e_2d_2\cos{\theta_2}}{1+e_2\cos{\theta_2}}}\right)} = \arctan{\left(\frac{-\frac{e_1d_1\sin{\theta_1}}{1+e_1\cos{\theta_1}}}{2j_1-\frac{e_1d_1\cos{\theta_1}}{1+e_1\cos{\theta_1}}}\right)} + 2\pi k$$
Simplifying... $$\tan{\left(\frac{2\pi t}{T_1}\right)}e_2d_2(4j_1-e_1d_1) = tan{\left(\frac{2\pi t}{T_2}\right)}e_1d_1(4j_2-e_2d_2)$$
At which point, I got stuck. WolframAlpha does not give me a clear answer for $t$, and if I substitute the example gif variables in this formula, WolframAlpha gives 5 possible $t$, of which 3 contain an imaginary number $i$. Which doesn't make sense since time doesn't have 2 dimensions.
Questions
- Am I doing something wrong?
- Is there a more elegant solution to this problem that I'm not aware of?
- Why is WolframAlpha giving me imaginary numbers?

With respect to ellipse center the angles using standard notation $(a,b,c,e)$ major, minor axes, semi inter focal length distance and eccentricity respectively are found as:
$$ (x,y)= \dfrac{a(1-e^2)}{1+e\cos \theta} (\cos \theta, \sin \theta)$$
$$ (x',y')= \dfrac{-a(1-e^2)}{1+e\cos \theta} (\cos \theta, \sin \theta)$$
$$ \omega= \theta + \tan^{-1}\dfrac{y}{c-x}+ \tan^{-1}\dfrac{y}{2c-x}; \tag 1 $$
$$ \theta_1 + \tan^{-1}\dfrac{y_1}{c_1-x_1}+ \tan^{-1}\dfrac{y_1}{2c_1-x_1}$$
$$=\theta_2 + \tan^{-1}\dfrac{y_2}{c_2-x_2}+ \tan^{-1}\dfrac{y_2}{2c_2-x_2};\tag 2 $$
Eccentric anomaly $E$ at center is connected to the $\theta's$ from Kepler's Equation:
$$ y_1/x_1=\tan E_1=\sqrt{\dfrac{1+e_1}{1-e_1}}\tan \theta_1;\;$$ $$ y_2/x_2=\tan E_2=\sqrt{\dfrac{1+e_2}{1-e_2}}\tan \theta_2 ;\tag 3$$
From Kepler's third law
$$ \frac{T_1^2}{T_2^2}= \frac {a_1^3}{a_2^3} \tag4 $$
Time dependence is from the classical Newton planetary orbits ( tangential direction & Coriolis accelerations and normal direction the gravity & centrifugal accelerations ) withgiven perigee start boundary conditions. we can calculate all time dependent dynamical quantities either using elliptic integrals in closed form or numerically. Finding roots of angles establishes times of clock angles.
Fuller answer could follow after OP's clarifying reply to my comments.