This situation is somewhat more complicated than you’ve assumed in your question. The center of the ellipse that’s the projection of the circle will in general not be the projection of the circle’s center. It will “drift” away from this point. This is because for a finite camera, i.e., one positioned a finite distance from the image, the projection is asymmetric: rays from the camera to the edge of the circle that’s tilted away from the camera make a smaller angle with the camera’s axis than do the corresponding rays on the opposite edge, so the foreshortening parallel to the $y$-axis is greater on that side. This asymmetry also means that uniformly stretching the magic paper doesn’t correspond to a uniform scaling of the image (glass) plane, and vice-versa. We’ll come back to this later.
Aside from the obvious dependency on the tilt angle $\theta$, for a camera positioned a finite distance $f\gt0$ from the image plane, the size, shape and location of the circle’s image also depend on the distance $d\gt0$ of the circle’s center from the image plane and radius $r$ of the circle. Note that even though we’re not given this radius explicitly, we can compute it from the two given points: $r=\sqrt{(x-x_0)^2+(y-y_0)^2}$. It’s reasonable to restrict $\theta$ to the open interval $\left(-\frac\pi2,\frac\pi2\right)$: When $\theta=\pm\frac\pi2$ we’re looking at the paper edge on, and beyond that we’re looking at the back of the paper.
It seems easiest to me to proceed in two stages: compute the raw projection of the circle and then adjust it to fit the constraints. We’ll do this by computing a homography—a planar perspective projection—between the two planes. W.l.o.g. assume that $(x_0,y_0)=(0,0)$ since we can always translate to make it so. Place the glass on the $x$-$y$ plane. Working in homogeneous coordinates, the camera is at $\mathbf V=[0,0,f,1]^T$ and the circle’s center at $[0,0,-d,1]^T$. A normal to the paper’s plane is $[0,\cos\theta,\sin\theta,0]^T$. The image plane is then $\mathbf\Pi_{\text{im}}=[0,0,-1,0]^T$ and the paper’s plane is $\mathbf\Pi_{\text{src}}=[0,\sin\theta,-\cos\theta,-d\cos\theta]^T$, although we won’t make direct use of the latter. The matrix $\mathcal P$ of the projection onto the image plane is $$\mathcal P=\mathbf V\mathbf\Pi_{\text{im}}^T-\mathbf V^T\mathbf\Pi_{\text{im}}\mathcal I_4=\begin{bmatrix}f&0&0&0\\0&f&0&0\\0&0&0&0\\0&0&-1&f\end{bmatrix}.$$ To complete the unadjusted homography we choose convenient bases for the two planes. For the image plane, we’ll simply drop the $z$-coordinate. For the source plane, we place the origin at the circle’s center, take the “$x$”-axis parallel to the world $x$-axis and use the normal for the “$z$”-axis. Putting these together with $\mathcal P$ produces the homography $$\mathcal H=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&0&1\end{bmatrix}\begin{bmatrix}f&0&0&0\\0&f&0&0\\0&0&0&0\\0&0&-1&f\end{bmatrix}\begin{bmatrix}1&0&0\\0&\cos\theta&0\\0&\sin\theta&-d\\0&0&1\end{bmatrix} = \begin{bmatrix}f&0&0\\0&f\cos\theta&0\\0&-\sin\theta&f+d\end{bmatrix}.$$
The general equation of a plane conic can be written in matrix form as $[x,y,1]\mathcal C[x,y,1]^T=0$, where $\mathcal C$ is a symmetric $3\times3$ matrix. This equation has the same form when “homogenized” by replacing $x$ and $y$ by $x/w$ and $y/w$, respectively, and we can identify conics with these matrices. They are covariant: if points transform according to $\mathbf p'=\mathcal M\mathbf p$, then $\mathcal C$ transforms as $\mathcal M^{-T}\mathcal C\mathcal M^{-1}$. Our circle is represented by the matrix $\mathcal C=\operatorname{diag}(1,1,-r^2)$ and its image is: $$\mathcal C'=\mathcal H^{-T}\mathcal C\mathcal H^{-1}=\begin{bmatrix} (f+d)^2&0&0 \\ 0&(f+d)^2\sec^2\theta-r^2\tan^2\theta&-r^2f\tan\theta \\ 0&-r^2f\tan\theta&-r^2f^2 \end{bmatrix}.$$ (I’ve multiplied through by $f^2(f+d)^2$, which doesn’t change the conic.) $\det\mathcal C'\lt0$, so this is an ellipse when $(f+d)^2\gt r^2\sin^2\theta$, an hyperbola when $(f+d)^2\lt r^2\sin^2\theta$, while equality produces a parabola. These conditions correspond to zero, two and one intersections, respectively, of the tilted circle with the plane $z=f$. This ellipse is axis-aligned, but its center is not at the origin. $\mathcal H$ preserves colinearity, so the center of this ellipse is the midpoint of the points $\mathcal H[0,\pm r,1]^T$, $$\mathbf c'=\left[0,{rf\sin\theta\over(f+d)^2-r^2\sin^2\theta}r\cos\theta,1\right]^T.$$
By linearity, a translation in one plane corresponds to a translation in the other plane, so the projected image’s center can be adjusted by “sliding” the paper a bit in the direction of its “$y$”-axis. On the other hand, because of the asymmetry in the projection noted at the beginning, uniformly stretching the paper to hit the target point will change the shape of the ellipse. Not only that, because the center of the projection depends nonlinearly on the circle’s size, this stretch will also affect the shift needed to place the center where it’s supposed to go. I don’t think that’s what you really had in mind when you described the “magic paper,” so I’ll go through how to make the necessary adjustments after projection.
Since we already have the ellipse center $\mathbf c'$, it’s easiest to first translate and then scale. The effect of this translation is easy to compute: the translated ellipse is in standard position, so the resulting matrix is diagonal, with the first two elements the same as those along the main diagonal of $\mathcal C'$, while the lower-right element becomes $\mathbf c'^T\mathcal C'\mathbf c'$. We can then normalize this matrix so that the lower-right element is $-1$, after which the semi-axis lengths can be read off as the reciprocal square roots of the other two diagonal elements. That is, $$\begin{align}a^2 &= -{\mathbf c'^T\mathcal C'\mathbf c' \over \mathcal C'_{1,1}} \\ b^2 &= -{\mathbf c'^T\mathcal C'\mathbf c' \over \mathcal C'_{2,2}}.\end{align}$$ However, before going to the trouble of calculating these values, remember that this ellipse will be scaled, so we really only care about its eccentricity, which depends on the ratio of $a^2$ to $b^2$. Thus, we can simply use the diagonal elements $\mathcal C'_{1,1}$ and $\mathcal C'_{2,2}$ directly without computing $\mathbf c'^T\mathcal C'\mathbf c'$ (or, for that matter, the center $\mathbf c'$ itself) and even discard other common factors to get $$a^{-2}=1 \\ b^{-2}=\sec^2\theta-{r^2\over(f+d)^2}\tan^2\theta.$$
(To satisfy my curiosity, I computed the actual values, anyway: $$a^2 = {f^2\over(f+d)^2-r^2\sin^2\theta}r^2 \\ b^2 = \left({f(f+d)\over(f+d)^2-r^2\sin^2\theta}\right)^2r^2\cos^2\theta.$$ As $f\to\infty$, i.e., as the projection becomes parallel, $a\to r$, $b\to r\cos\theta$ and $\mathbf c'\to[0,0,1]^T$.)
To find the necessary scale factor so that this translated ellipse passes through $\mathbf p'=[x_p,y_p,1]^T$, observe that scaling the ellipse ${x^2\over a^2}+{y^2\over b^2}=1$ by a factor of $s>0$ changes the equation to ${x^2\over a^2}+{y^2\over b^2}=s^2$, therefore the correct scale factor is $$s=\sqrt{{x_p^2\over a^2}+{y_p^2\over b^2}}$$ and the equation of the adjusted ellipse is $${x^2\over a^2}+{y^2\over b^2}={x_p^2\over a^2}+{y_p^2\over b^2}$$ with $a^2$ and $b^2$ as above. Substituting for $a^2$ and $b^2$ and rearranging finally gives the sought-after equation: $$(x^2-x_p^2)+(y^2-y_p^2)\left(\sec^2\theta-{r^2\over(f+d)^2}\tan^2\theta\right)=0.\tag{*}$$ In the limit as $f\to\infty$, this becomes $(x^2-x_p^2)+(y^2-y_p^2)\sec^2\theta$. Note that the final adjusted ellipse doesn’t so much depend on the actual values of $r$, $f$ and $d$ as on the ratio of the size of the circle to the distance of the camera from its center. If this distance is large relative to the radius of the circle, or the tilt angle very large or small, the limit ellipse is a good approximation to the true image.
In matrix terms, the final ellipse is obtained via the transformation matrix $$\begin{bmatrix}s&0&-sc'_x\\0&s&-sc'_y\\0&0&1\end{bmatrix}\mathcal H = \begin{bmatrix}sf&0&0 \\ 0&s(f\cos\theta+c'_y\sin\theta)&-sc'_y(f+d) \\ 0&-\sin\theta&f+d \end{bmatrix}.$$ Here, $s$ should be computed using the true values of $a$ and $b$. This matrix allows you to map the entire source plane, not just the circle, if so desired.
If you want to do the scaling on the paper side instead, that will entail a similar computation. Instead of starting out with $\mathcal H$, however, you’ll use $$\mathcal H'=\begin{bmatrix}sf&0&0\\0&sf\cos\theta&0\\0&-s\sin\theta&f+d\end{bmatrix}.$$ The center of the ellipse and its semi-axis lengths are found as above, but they will all depend on the unknown scale factor $s$. Plugging the coordinates of $\mathbf p'$ into the resulting ellipse equation will let you solve for this scale factor.
