I am trying to compute a part of the projected area in the x-y plane (base area) of an ellipsoidal cap rotated in the x-y plane by an angle $\beta$. So I am taking a point on the rotated ellipsoid given by its parametric form:
$$x = (a \cos(\theta) \cos(\beta)-b \sin(\theta) \sin(\beta)) \sin(\phi)$$ $$y = (a \cos(\theta) \sin(\beta)+b \sin(\theta) \cos(\beta)) \sin(\phi)$$ $$z = c \cos(\phi)$$
where the axis size in x coordinate is 'a', the axis size in y-coordinate is 'b' and the axis size in z-coordinate is 'c'. $(\theta)$ is the angle with x axis in x-y plane and $(\phi)$ is the angle with z axis in the spherical coordinates:
Now for full projected area in the x-y plane I double integrate the following from $0$ to $2\pi$ and $0$ to $\phi_0$, where $\phi_0$ is the angle with z-axis at the cap height:
$$A_{xy} = \iint \mid{dA_{xy}}\mid = \iint dx \times dy = \iint f(\theta,\phi) d\theta d\phi $$
I integrate the area over an absolute value of $dA_{xy}$, otherwise I get $0$. Following the steps I differentiate x and y partially over $\theta$ and $\phi$
$$dx = -(a \sin(\theta) \cos(\beta)+b \cos(\theta) \sin(\beta)) \sin(\phi) d\theta + (a \cos(\theta) \cos(\beta)-b \sin(\theta) \sin(\beta)) \cos(\phi) d\phi$$
$$dy = -(a \sin(\theta) \sin(\beta)-b\cos(\theta)\cos(\beta)) \sin(\phi))d\theta +(a \cos(\theta) \sin(\beta)+b\sin(\theta)\cos(\beta))\cos(\phi) d\phi$$
$$dz = -c\sin(\phi)d\phi$$
$$ dA_{xy} = dx*dy = \frac{1}{2}sin(2\phi)(ab\cos(2\theta)cos(2\beta)-(\frac{a^2+b^2}{2}sin(2\theta)sin(2\beta))d\theta d\phi$$
After integration my total projected area in the x-y plane changes with the angle $(\beta)$. It should stay constant irrespective of the angle of rotation.Is my approach correct? Or is there an error in my method which I cannot find? Please help.
You were doing fine up through your calculation of $dx$ and $dy$, but then took a wrong turn when you computed the resulting volume element. I’m not really sure what you did to end up with the expression that you have for what you call $dA_{xy}$, but you can’t just multiply the expressions that you have for $dx$ and $dy$ to compute it. You’re essentially performing a coordinate transformation on $dA$, so you have to use the change-of-coordinate formula $$dx\,dy = \det{\partial(x,y) \over \partial(\phi,\theta)}d\theta\,d\phi,$$ with the Jacobian determinant as the conversion factor. If you work through this correctly, you’ll end up with a somewhat messy expression that involves neither $\theta$ nor $\beta$.
However, there’s a simpler way to perform this area calculation. If you fix $\phi$ and set $z$ to zero, you have a parameterization of the boundary of the projection. This suggests using Green’s theorem to turn $\iint_RdA$ into a single integral over the boundary. There are many antiderivatives of $1$, but experience suggests that the most convenient one for this problem is $\frac12(x\,dy-y\,dx)$. This is the area of the “infinitesimal” triangle with sides $(x,y)$ and $(dx,dy)$ that approximates the area swept out by a radial vector. In this case, the integrand eventually simplifies to $\frac12ab\sin^2\phi_0$, which results in an area of $\pi ab\sin^2\phi_0$ as one would expect.