Option 1:
The projection of the intersection of $z=x^2+y^2$ and $z=y$ in the $xy$ plane is $x^2+y^2=y$, i.e., the circle centered at $(0,1/2)$ with radius $1/2$, or in polar coordinates, $r=\sin \theta$. So you can parametrize the surface $S$ as follows:
$$
\begin{cases}
x=x \\
y=y \quad \quad \quad \quad \mbox{with }(x,y)\in D=\{(r,\theta)\;|\; 0 \le \theta \le \pi, \; 0 \le r \le \sin \theta \}\\
z=x^2+y^2
\end{cases}
$$
It follows that
$$
\Phi = \iint_S \nabla \times \vec{v} \; dS = \iint_D \nabla \times \vec{v}(x,y)\cdot \pmatrix{1\\0\\2x}\times \pmatrix{0\\1\\2y}\; dA = \iint_D \pmatrix{1\\2\\1} \cdot \pmatrix{-2x\\-2y\\1}\; dA \\= \iint_D 1-2x-4y \;dA
$$
This matches your answer, except for the $4$ coefficient for $y$. If I am not mistaken, $ \nabla \times \vec{v}=(1,2,1)$. Now what you need to do is integrate over $D$ (or $A$ with your notation), as defined above.
Switching to polar coordinates yields
$$
\Phi = \int_0^{\pi}\int_0^{\sin \theta}(1-2r \cos \theta-4r\sin\theta) r dr d\theta = -\frac{\pi}{4}
$$
Option 2:
Alternatively, if you are familiar with the divergence theorem:
$$
\Phi = \iint_S \nabla \times \vec{v} \; dS = \underbrace{\iiint_V \nabla\cdot \nabla \times \vec{v} \; dV}_{=0} - \iint_{S_2} \nabla \times \vec{v} \; dS
$$
where $S_2$ is the surface that closes $S$, i.e., the part of the plane $z=y$ on top of $D$, which we can parametrize with
$$
\begin{cases}
x=x \\
y=y \quad \quad \quad \quad \mbox{with }(x,y)\in D=\{(r,\theta)\;|\; 0 \le \theta \le \pi, \; 0 \le r \le \sin \theta \}\\
z=y
\end{cases}
$$
It follows that
$$
\Phi = - \iint_{S_2} \nabla \times \vec{v} \; dS = - \iint_D \nabla \times \vec{v}(x,y)\cdot \pmatrix{1\\0\\0}\times \pmatrix{0\\1\\1}\; dA = - \iint_D \pmatrix{1\\2\\1} \cdot \pmatrix{0\\-1\\1}\; dA = \iint_D dA = \frac{\pi}{4}
$$
When we use this theorem, we conventionally parametrize the surface "outwards", i.e., in the opposition direction as in Option $1$, hence the sign difference.
Option 3:
Last but not least, you can use Stoke's theorem:
$$
\Phi = \iint_S \nabla \times \vec{v} \; dS = \oint_C \vec{v} \cdot d\vec{r}
$$
where $C$ is the curve at the intersection between the paraboloid and the plane. Using the parametrization proposed by @levap in the question that you asked earlier :
$$
\begin{cases}
x=\frac{\cos t}{2} \\
y=\frac{1+\sin t}{2} \quad \quad \quad \quad \mbox{with } t\in [0,2\pi] \\
z=\frac{1+\sin t}{2}
\end{cases}
$$
It follows that
$$
\Phi = \oint_C \vec{v} \cdot d\vec{r} = \int_0^{2\pi} \frac{1}{4} \pmatrix{2(1+\sin t)\\ 1+ \cos t \\ 1 + \sin t}\cdot \pmatrix{-\sin t\\ \cos t \\ \cos t}\; dt = -\frac{\pi}{4}
$$
Isn't Calculus fun ? :)
For any point $(x, y, z)$ on the green plane, you have
$$
\tan \theta = \frac{z}{x}.
$$
hence $z = x \tan(\theta)$.
Letting $t = \tan \theta$, to save a little typing, this means that your third line of code should be something like
function z_cylinder_plane_intersection(y_angle, z_angle, radius) =
(radius * cos(z_angle)) * tan(theta)
Best Answer
If I've understood well you need to find the equation of the intersection curve between a torus and a cylinder (I'll call it $tc$) defined by the following equations:
Torus: (1) $\,\,\,(x^2+y^2+z^2+R^2-r^2)^2=4 R^2 (x^2+y^2)$
Cylinder: (2) $\,\,\,z^2 + (y - R)^2= r^2 $
(BTW, according to your figure are you sure it's not $z^2 + (y + R)^2= r^2 $ ? anyway little changes in what follows)
You can isolate $z^2$ in the equation (2) and put it in the equation (1).
With some simplification you'll get
$x^2(x^2+4Ry-4R^2)=0$
That is
(3) $y = - \frac{{{x^2}}}{{4R}} + R$
This means that the projection of the $tc$ curve on the $x$y plane is a simple parabola.
After that you can put equation (3) into equation (2) and get (4): ${z^2} = - \frac{{{x^4}}}{{16{R^2}}} + {r^2}$
That's the projection of the $tc$ curve in the $xz$ plane. It's a quartic curve.
Your $tc$ curve in space can then be defined as the intersection of the surfaces defined by equations (3) and (4) where the first is the surface generated by translating the parabola of the $xy$ plane along the $z$ axis and the second is generated by translating the quartic curve of the $xz$ plane along the $y$ axis.
The parametric equations of the $tc$ curve are then
$\begin{array}{*{20}{l}} {x = t}\\ {y = - \frac{{{t^2}}}{{4R}} + R}\\ {{z} = \pm \sqrt { - \frac{{{x^4}}}{{16{R^2}}} + {r^2}} } \end{array}$
If you're interested you can also derive the projection of the intersection curve in the $yz$ plane. I think you'll get a semicircle.
I hope this helps
[I suggest to add the tags "solid-geometry" and "3D" to your original post]
Luca